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Abstract 


This  work  develops  an  adaptive  concurrent  multi-level  computational  model  for  multi-scale  analysis  of 
composite  structures  undergoing  damage  initiation  and  growth  due  to  microstructural  damage  induced  by 
debonding  at  the  fiber-matrix  interface.  The  model  combines  macroscopic  computations  using  a  continuum 
damage  model  developed  in  a  preceding  paper  [84]  with  explicit  micromechanical  computations  of  stresses 
and  strain,  including  explicit  debonding  at  the  fiber-matrix  interface.  The  macroscopic  computations  are 
done  by  conventional  FEM  models  while  the  Voronoi  cell  FEM  is  used  for  micromechanical  analysis.  Three 
hierarchical  levels  of  different  resolution  adaptively  evolve  in  this  to  improve  the  accuracy  of  solutions  by 
reducing  modeling  and  discretization  errors.  They  levels  include:  (a)  level-0  of  pure  macroscopic  analy¬ 
sis  using  a  continuum  damage  mechanics  (CDM)  model;  (b)  level- 1  of  asymptotic  homogenization  based 
macroscopic-microscopic  RVE  modeling  to  monitor  the  breakdown  of  continuum  laws  and  signal  the  need 
for  microscopic  analyses;  and  (c)  level-2  regions  of  pure  micromechanical  modeling  with  explicit  depiction  of 
the  local  microstructure.  Two  numerical  examples  are  solved  to  demonstrate  the  effectiveness  and  accuracy 
of  the  multi-scale  model.  A  double  lap  bonded  composite  joint  is  modeled  for  demonstrating  the  model’s 
capability  in  handling  large  structural  problems. 

For  micromechanical  analysis,  the  traditional  Voronoi  cell  finite  element  model(VCFEM)  has  been  im¬ 
proved  for  studying  the  interfacial  debonding.  A  new  extended  Voronoi  cell  finite  element  model(X-VCFEM) 
has  been  developed  for  modeling  interfacial  debonding  with  arbitrary  matrix  cohesive  cracking  in  fiber- 
reinforced  composites.  To  describe  the  onset  and  growth  of  damage  along  the  fiber-matrix  interface,  normal 
and  tangential  cohesive  zone  models  are  coupled  into  VCFEM.  It  is  shown  that  the  initiation  and  espe¬ 
cially  propagation  of  debonding  depends  not  only  on  the  total  cohesive  energy,  but  also  on  the  shape  of  the 
traction-displacement  curve.  The  model  is  also  used  to  study  the  influence  of  various  local  morphological 
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parameters  on  damage  evolution  by  interfacial  debonding.  A  special  function  of  various  geometric  parame¬ 
ters  is  developed  to  predict  the  location  of  debonding  in  microstructures  with  varying  morphology. 

Next  an  extended  Voronoi  cell  finite  element  model  (X-VCFEM)  is  developed  for  modeling  multiple 
cohesive  crack  propagation  in  brittle  materials.  The  cracks  are  modeled  by  a  cohesive  zone  model  and  their 
incremental  directions  and  growth  lengths  are  determined  in  terms  of  the  cohesive  energy  near  the  crack  tip. 
Extension  to  VCFEM  is  achieved  through  enhancements  in  stress  functions  in  the  assumed  stress  hybrid 
formulation.  In  addition  to  polynomial  terms,  the  stress  functions  include  branch  functions  in  conjunction 
with  level  set  methods,  and  multi-resolution  wavelet  functions  in  the  vicinity  of  crack  tips.  Comparison  of 
X-VCFEM  simulation  results  with  results  in  literature  for  several  fracture  mechanics  problems  validates  the 
effectiveness  of  X-VCFEM.  Effect  of  stereographic  features  such  as  size  and  distribution  of  heterogeneities 
on  damage  evolution  in  random  microstructures  are  also  discussed. 

Finally,  in  order  to  study  the  interaction  between  interface  debonding  and  cohesive  matrix  cracking,  a 
criterion  based  on  cohesive  models  is  proposed  to  assess  the  crack  penetrating  into  matrix  from  the  interface 
and  is  validated  by  numerical  examples. 
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Chapter  2 


Introduction 


The  commercial  use  of  advanced  fiber  reinforced  composites  in  structural  components  has  considerably 
increased  due  to  their  superior  thermo-mechanical  properties.  Based  on  design  requirements,  they  are 
engineered  to  yield  high  strength  or  stiffness  to  weight  ratios,  resulting  in  a  tremendous  advantage  over 
conventional  materials.  Despite  mechanical  property  enhancements,  the  presence  of  fibers  in  composites 
often  has  adverse  effects  on  their  failure  properties.  Both  debonding  at  the  fiber-matrix  interface  and  crack 
in  matrix  are  major  micromechanical  damage  phenomena,  responsible  for  deterring  the  overall  properties 
and  resulting  in  diminished  structural  integrity.  The  search  for  a  rational  design  procedure  to  select  optimal 
composite  microstructures  provides  a  compelling  reason  for  accelerated  development  of  methods  relating  the 
microstructure  to  the  material’s  mechanical  behavior  and  failure  characteristics. 

Analysis  of  composite  materials  with  microstructural  heterogeneities  is  conventionally  done  with  macro¬ 
scopic  properties  obtained  by  homogenizing  response  functions  in  the  representative  volume  element  (RVE) 
from  microscopic  analyses  at  smaller  length  scales.  While  these  “bottom-up”  homogenization  models  are 
efficient  and  can  reasonably  predict  macroscopic  or  averaged  behavior,  such  as  stiffness  or  strength,  they 
have  limited  predictive  capabilities  with  problems  involving  localization,  failure  or  instability.  Assumptions 
of  macroscopic  uniformity  and  RVE  periodicity,  the  two  basic  requirements  of  homogenization,  break  down 
under  these  circumstances.  The  uniformity  assumption  ceases  to  hold  in  critical  regions  of  high  local  solution 
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gradients,  such  as  near  free  edges,  interfaces,  material  discontinuities  or  evolving  damage.  RVE  periodicity, 
on  the  other  hand,  is  unrealistic  for  non-uniform  microstructures,  e.g.  in  the  presence  of  clustering  of  hetero¬ 
geneities  or  microscopic  damage.  Even  with  a  uniform  phase  distribution  in  the  microstructure,  the  evolution 
of  localized  stresses,  strains  or  damage  path  can  violate  the  periodicity  conditions.  Problems  like  this  have 
been  effectively  tackled  by  multi-scale  modeling  methods  e.g.  in  [81,  33,  50,  75,  74,  90,  89,  83,  82, 101,  120,  99]. 
Multi-scale  analyses  methods  can  be  broadly  classified  into  two  classes.  The  first  is  known  as  ’’hierarchical 
models”  [33,  50,  101,  99]  in  which  information  is  passed  from  lower  to  higher  scales,  usually  in  the  form 
of  material  properties.  The  hierarchical  homogenization  models  assume  periodic  representative  volume  ele¬ 
ments  (RVE)  in  the  microstructure  and  uniformity  of  macroscopic  field  variables.  The  second  class,  known 
as  “concurrent  methods”  [89,  75,  74,  90,  83,  82,  120],  implement  sub-structuring  and  simultaneously  solve 
different  models  at  regions  with  different  resolutions  or  scales. 

The  two-way  coupling  of  scales  enabled  in  the  concurrent  methods  is  suitable  for  problems  involving 
localization,  damage  and  failure.  Macroscopic  analysis,  using  bottom-up  homogenization  in  regions  of  rel¬ 
atively  benign  deformation,  enhances  the  efficiency  of  the  computational  analysis.  As  a  matter  of  fact,  it 
would  be  impossible  to  analyze  large  structural  regions  without  the  advantage  of  a  continuum  model  based 
macroscopic  analysis.  On  the  other  hand,  the  top-down  localization  process  cascading  down  to  the  mi¬ 
crostructure  in  critical  regions  of  localized  damage  or  instability  for  pure  microscopic  analysis,  is  necessary 
for  accurately  predicting  the  damage  path.  These  microscopic  computations,  depicting  the  real  microstruc¬ 
ture  are  often  complex  and  computationally  prohibitive.  Hence,  a  concurrent  setting  makes  such  analyses 
feasible,  provided  the  ”zoom-in”  regions  are  kept  to  a  minimum.  The  adaptive  multi-level  models,  promoted 
in  [75,  74,  90,  83,  82,  120],  are  attempts  to  achieve  this  objective,  with  the  adaptivity  motivated  from  phys¬ 
ical  and  mathematical  perspectives.  However,  there  is  a  paucity  of  such  studies  in  the  literature  involving 
material  nonlinearity  and  evolving  microstructural  damage.  In  their  previous  studies,  Ghosh  and  coworkers 
have  proposed  adaptive  multi-level  analysis  using  the  microstructural  Voronoi  cell  FEM  model  for  modeling 
elastic-plastic  composites  with  particle  cracking  and  porosities  in  [89],  and  for  elastic  composites  with  free 
edges  and  stress  singularities  in  [83,  82]. 

In  this  work,  we  have  derived  and  computationally  modeled  an  anisotropic  continuum  damage  mechanics 
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(CDM)  model  for  unidirectional  fiber-reinforced  composites  undergoing  interfacial  debonding  from  by  using 
homogenization  theory.  The  CDM  model  homogenizes  the  damage  incurred  through  initiation  and  growth 
of  interfacial  debonding  in  a  microstructural  RVE  with  nonuniform  distribution  of  fibers.  Additionally, 
arbitrary  loading  conditions  are  also  effectively  handled  by  this  model.  The  CDM  model  is  then  used  in  an 
adaptive  concurrent  multi-level  computational  model  to  analyze  multi-scale  evolution  of  damage.  Damage 
by  fiber-matrix  interface  debonding,  is  explicitly  modeled  over  extended  microstructural  regions  at  critical 
locations  [37,  58].  The  adaptive  model  addresses  issues  of  efficiency  and  accuracy  through  considerations  of 
physically-based  modeling  errors. 

The  adaptive  multi-level  model  consists  of  three  levels  of  hierarchy  viz.  level-0,  level-1  and  level-2),  which 
evolve  in  sequence.  The  continuum  damage  model  developed  in  [84]  is  used  for  level-0  computations.  The 
level-1  domain  is  used  as  a  ‘swing  region’  to  establish  criteria  for  switching  from  macroscopic  to  microscopic 
calculations.  Physical  criteria  involving  variables  at  the  macroscopic  and  microstructural  RVE  levels,  trigger 
switching  from  pure  macroscopic  to  pure  microscopic  calculations,  i.e.  the  level  —  0  — >  level  —  1  — *  level  —  2. 
A  transition  layer  is  placed  between  the  level  —  1  and  microscopic  level  —  2  domains  for  smooth  transition 
from  one  scale  to  the  next. 

All  computations  in  the  composite  microstructure  with  explicit  representations  of  the  fiber  and  matrix 
phases  are  done  with  the  Voronoi  cell  finite  element  model  or  VCFEM  [37,  58].  Accurate  interface  debond¬ 
ing  analysis  is  a  difficult  task  due  to  the  fact  that  morphological  and  constitutive  complexities  govern  its 
initiation  and  growth.  The  present  study  is  aimed  at  understanding  the  effect  of  microstructure  on  interface 
decohesion  induced  damage  evolution  in  multi-fiber  microstructures.  The  effect  of  debonded  interface  on  the 
mechanical  properties  have  been  studied  by  several  investigators  e.g.  [13,  44,  77],  using  simplified  models  for 
representing  imperfect  conditions  through  traction  discontinuities.  The  propagation  of  interface  cracking  has 
been  successfully  modeled  by  a  number  of  researchers  using  the  cohesive  volumetric  finite  element  methods. 
Among  the  important  contributions  to  the  field  of  damage  evolution  by  normal  and  tangential  separation 
are  those  by  Needleman  [71,  72,  73],  Tvergaard  [105,  104]  ,  Allen  et.  al.  [4,  62],  Lissenden  et.  al.  [61], 
Geubelle  et.  al.  [35],  Walter  et.  al.  [109]  among  others.  Yuan  et.al.  [115]  and  Achenbach  and  Zhu  [1]  have 


10 


evaluated  microscopic  and  macroscopic  responses  for  imperfect  and  debonded  interfaces  using  finite  element 
models.  A  majority  of  these  studies  have  used  unit  cell  models  with  periodic  repetition  of  single  cells.  These 
models  provide  valuable  insights  into  the  microstructural  damage  processes.  However  many  of  these  are 
ineffective  in  predicting  the  interaction  between  fibers,  effects  of  clustering,  alignment,  relative  sizes  etc., 
that  are  critical  to  the  failure  process  in  the  microstructure.  Zhong  and  Knauss  [117]  have  proposed  a  hybrid 
discrete-continuum  approach  in  which  discrete  particle  interactions  with  damage  evolution  are  modeled,  ac¬ 
counting  for  particle  size  and  spacing. 

Another  important  damage  phenomenon  in  composites  is  crack  propagation  in  brittle  matrix.  The  dif¬ 
ficulty  in  obtaining  analytical  solutions  for  many  of  these  problems,  especially  those  associated  with  crack 
propagation,  has  prompted  the  use  of  numerical  methods  like  the  finite  element  method  for  the  determina¬ 
tion  of  fracture  mechanics  parameters  such  as  stress  intensity  factors,  energy  release  rates  and  J-integrals, 
crack  tip  stresses  and  opening  displacements.  Analyses  using  conventional  finite  element  method  require 
very  high  density  mesh  to  overcome  limitations  of  severe  pathological  mesh  dependence  near  the  crack  tips. 
Convergence  is  very  slow  since  stress  singularities  are  not  accounted  for  in  the  element  formulation.  For 
improving  the  computational  efficiency  through  better  representation  of  the  crack  tip  singularity,  a  number 
of  different  methods  have  been  proposed.  The  superposition  method  [113,  112]  has  introduced  the  superpo¬ 
sition  of  singular  terms  to  the  finite  element  interpolations.  The  singular  element  method  introducing  the 
quarter-point  elements  [7,  8,  47,  48]  near  the  crack  tip  has  been  developed  to  yield  reasonably  accurate  crack 
tip  parameters  and  displacements.  As  an  alternative  to  the  displacement  based  finite  element  models,  hybrid 
singular  elements  have  been  proposed  in  [103,  102,  60,  78,  53,  116].  Also  termed  as  super-elements,  these 
elements  accommodate  crack  tip  singularity  through  interpolation  functions  that  account  for  stress  intensity 
factors  using  classical  elasticity  theory.  Furthermore,  numerical  analysis  and  simulation  of  the  growth  of 
multiple  cracks  in  materials  is  a  challenging  enterprise  due  to  morphological  and  constitutive  complexities 
that  govern  its  growth.  Even  a  very  high  density  mesh  cannot  overcome  pathological  mesh  dependence 
near  the  crack  tips  and  avoid  biasing  the  direction  of  crack  propagation.  The  difficulties  aggravate  in  the 
presence  of  multiple  cracks,  due  to  their  interaction  with  each  other.  Various  methods  have  been  proposed 
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for  improving  the  effectiveness  of  computational  methods  in  modeling  cracks.  While  most  of  these  analyses 
are  limited  to  stationary  cracks,  it  is  only  recently  that  effective  methods  of  analysis  of  crack  propagation 
are  being  proposed. 

With  increasing  power  of  computational  modeling  and  hardware,  the  cohesive  zone  models  [71,  72,  73, 
104,  34,  35,  42,  76]  have  emerged  as  important  tools  for  modeling  crack  propagation  in  homogeneous  and 
heterogeneous  materials.  In  these  models,  interfaces  of  similar  and  dissimilar  materials  are  treated  as  zero 
thickness  non-linear  springs.  Interfacial  traction  is  specified  as  nonlinear  functions  of  tangential  and  normal 
separations  across  the  interface  to  manifest  crack  evolution.  These  models  have  been  used  to  simulate  crack 
growth  between  elements  in  [16,  110,  42],  by  lacing  the  interface  between  contiguous  elements  with  cohesive 
springs.  The  use  of  a  highly  refined  computational  mesh,  especially  near  the  crack  tip  is  still  a  requirement, 
even  though  the  effect  is  mitigated  due  to  the  finite  crack  tip  stress  with  this  model.  Alternatively,  intra- 
element  enrichment  approaches,  based  on  the  incorporation  of  embedded  discontinuities  in  displacement  or 
strain  fields  have  been  proposed  ([52]),  which  eliminates  mesh  dependent  prediction  of  the  evolving  crack 
path,  and  hence  the  need  for  remeshing.  The  extended  FEM  or  X-FEM  [10,  9,  11,  12,  28,  66,  67]  is  a  powerful 
recent  addition  to  this  family  of  intra-element  enrichment.  Cohesive  crack  propagation  has  been  modeled  in 
this  work  by  using  the  partition  of  unity  concept  to  incorporate  local  enrichment  functions  that  allows  the 
preservation  of  the  general  displacement  based  FEM  formalism. 

Stress-based  finite  element  methods  have  had  considerable  success  when  stress  fields  are  of  interest  in 
the  analysis  [103,  102].  Within  this  general  formalism,  the  Voronoi  cell  finite  element  method  (VCFEM) 
has  been  developed  in  [37,  68,  88,  38,  87,  86,  58]  for  micromechanical  analysis  of  arbitrary  heterogeneous 
microstructures.  The  method  can  effectively  overcome  requirements  of  large  degrees  of  freedom  in  conven¬ 
tional  finite  element  models.  Morphological  arbitrariness  in  dispersions,  shapes  and  sizes  of  heterogeneities, 
as  seen  in  real  micrographs  are  readily  modeled  by  this  method.  The  VCFE  model  naturally  evolves  by 
tessellation  of  the  microstructure  into  a  network  of  multi-sided  Voronoi  polygons.  Each  Voronoi  cell  with 
embedded  heterogeneities  (particle,  fiber,  void,  crack  etc.)  represents  the  region  of  contiguity  for  the  hetero- 
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geneity,  and  is  treated  as  an  element  in  VCFEM.  VCFEM  elements  are  considerably  larger  than  conventional 
FEM  elements  and  incorporate  a  special  assumed  stress  hybrid  FEM  formulation.  Incorporation  of  known 
functional  forms  from  analytical  micromechanics  substantially  enhances  its  convergence.  A  high  level  of 
accuracy  with  significantly  reduced  degrees  of  freedom  has  been  achieved  with  VCFEM.  Computational 
efficiency  is  therefore  substantially  enhanced  compared  to  conventional  displacement-based  FE  models.  Suc¬ 
cessful  applications  of  2D-small  deformation  VCFEM  have  been  made  in  thermo-elastic-plastic  problems 
of  composite  and  porous  materials  [68,  88],  An  adaptive  VCFEM  has  been  developed  in  [88],  where  opti¬ 
mal  improvement  is  achieved  by  h-p  adaptation  of  the  displacement  field  and  p-enrichment  of  the  stress  field. 

The  cohesive  crack  propagation  model  has  been  incorporated  in  VCFEM  in  [37,  58]  to  model  interface 
debonding  in  fiber  reinforced  composites.  However,  in  these  models,  the  debonding  or  crack  evolution  path 
is  along  the  interface  and  hence  the  cohesive  zone  regions  are  known  a-priori.  In  the  event  that  the  crack 
branches  off  into  the  matrix,  the  path  is  no  longer  pre-assessed  and  needs  to  be  determined  at  each  load 
increment,  consistent  with  the  local  state  of  stresses,  strains  and  morphology.  This  task  is  considerably  more 
challenging  since  a  slight  deviation  can  lead  to  completely  wrong  prediction. 

The  motivation  of  this  work  is  derived  from  the  need  to  create  a  robust  finite  element  method,  extended 
Voronoi  cell  finite  element  model  (X-VCFEM),  for  modeling  interface  debonding  with  arbitrary  crack  prop¬ 
agation  in  heterogeneous  materials.  This  is  an  essential  step,  prior  to  simulating  the  entire  microstructural 
failure  problem.  X-VCFEM  incorporates:  (a)  stress  discontinuities  across  the  cohesive  crack  through  branch 
functions  in  conjunction  with  level  set  methods,  (b)  crack  tip  stress  concentration  through  the  introduction 
of  multi-resolution  wavelet  functions  [41,  51,  80]  in  the  vicinity  of  the  crack  tip,  and  (c)  incremental  crack 
propagation  using  a  cohesive  energy  based  criterion  for  estimating  the  direction  and  length  of  the  incremental 
crack  advance. 
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2.1  Organization  of  this  Report 


The  report  is  divided  into  five  subsequent  chapters.  In  chapter  3,  variational  formulation  and  various 
aspects  of  the  computational  scheme  for  interface  debonding  in  composites  are  presented.  In  this  chapter, 
the  VCFE  model  is  also  used  to  study  the  influence  of  various  local  morphological  parameters  on  damage 
evolution  by  interface  debonding  and  the  sensitivity  of  the  debonding  process  to  spatial  distribution  and 
fiber  size.  Extensions  of  the  VCFEM  (X-VCFEM)  for  cohesive  crack  propagation  are  developed  in  Chapter 
4.  Numerical  validation  of  X-VCFEM  for  matrix  cracking  is  also  presented.  Based  on  the  preparation  of 
previous  two  chapters,  the  X-VCFEM  for  modeling  interface  debonding  with  matrix  cohesive  cracking  are 
developed  in  chapter  5.  Finally,  an  account  of  multi-scale  modeling  is  presented  with  critical  examples  The 
report  ends  with  a  discussion  and  conclusion  of  the  overall  effort  in  chapter  7. 
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Chapter  3 

Voronoi  cell  Finite  Element  Method 
for  Interfacial  Debonding  in 
Composite  Microstructures  with 
Morphological  Variations 

3.1  Introduction 

Debonding  at  the  fiber-matrix  interface  is  a  major  micromechanical  damage  phenomenon  in  composites, 
responsible  for  deterring  the  overall  properties  and  resulting  in  diminished  structural  integrity.  The  search 
for  a  rational  design  procedure  to  select  optimal  composite  microstructures  provides  a  compelling  reason 
for  accelerated  development  of  methods  relating  the  microstructure  to  the  material’s  mechanical  behavior 
and  failure  characteristics.  The  present  study  is  aimed  at  understanding  the  effect  of  microstructure  on 
interfacial  decohesion  induced  damage  evolution  in  multi-fiber  microstructures.  However  many  of  these  are 
ineffective  in  predicting  the  interaction  between  fibers,  effects  of  clustering,  alignment,  relative  sizes  etc., 
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with  VCFEM,  the  inclusion-matrix  interface  in  each  element  is  lined  with  a  series  of  cohesive  zone  springs, 
for  which  debonding  is  a  phenomenon  of  progressive  separation  across  an  extended  crack  tip  that  is  resisted 
by  cohesive  tractions  depicted  in  figure  3.1(b). 


3.2.1  Assumed  Stress  VCFEM  Formulation 

Each  Voronoi  cell  element  is  composed  of  the  matrix  phase  (flm),  the  inclusion  phase(Oc)  and  a  zero  thickness 
interface  region(ft{n  :  volume( fljn)  — *  0),  such  that  =  Qm  (J  |Jflin.  The  element  outer  boundary 
consists  of  the  prescribed  displacement  boundary  (rum),  prescribed  traction  boundary  (Ffm)  and  the  inter¬ 
element  boundary  (rm),  so  i.e.  dCle  =  rurn  (J  F(m  (J  Tm.  Compatible  displacement  conditions  apply  on 
dCle.  On  the  other  hand,  the  inner  matrix-inclusion  interface  (dQ.™ / dfl^)  in  each  element  is  facilitated  with 
incompatible  displacements  across  it  as  seen  in  figure  3.1(b).  In  order  to  describe  debonding  with  progressing 
deformation  through  decohesion,  the  interface  is  lined  with  a  set  of  node-pairs  with  nodes  belonging  to  the 
matrix  interface  ( dfl ™)  and  inclusion  interface  (dClf)  respectively.  d£lcc  has  an  outward  normal  nc  (=nm), 
while  ne  is  the  outward  normal  to  dfle.  In  the  incremental  assumed  stress  hybrid  VCFEM  formulation 
(37,  68,  88],  the  complementary  energy  functional  for  each  element  is  expressed  in  terms  of  increments  of 
stress  and  displacement  fields  as: 


ne(<7,  Aa,  u,  Au)  =  -  [  AB(am,Aam)dn  -  [  AB{ac,Aac)dn 
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Here  B  is  the  complementary  energy  density  and  the  superscripts  m  and  c  correspond  to  variables  associated 
with  the  matrix  and  inclusion  phases.  crm  and  <jc  are  the  equilibrated  stress  fields,  em  and  ec  the  corre¬ 
sponding  strain  fields  in  different  phases  of  each  Voronoi  element.  Also,  u,  um  and  uc  are  the  kinematically 
admissible  displacement  fields  on  d£lt,  dQ™  and  dflcc  respectively.  The  prefix  A  corresponds  to  increments. 
The  last  two  terms  provide  the  work  done  by  the  interfacial  tractions  Tm  =  T™ nm  +Ttm tm  due  to  interfacial 
separation  (um  —  uc),  where  T™  and  T(m  are  the  normal  and  tangential  components  that  are  described  by 
cohesive  laws.  The  total  energy  for  the  entire  composite  domain  is  obtained  by  adding  the  energy  functionals 
for  N  elements  as 

N 

n  =  £ne  (3.2) 

e=l 

Independent  assumptions  on  stress  increments  Act  are  made  in  each  of  the  element  phases  to  accommodate 
stress  jumps  across  the  interface.  A  convenient  method  of  deriving  equilibrated  stress  increments  in  each 
phase  is  through  the  introduction  of  stress  functions  &(x,  y),  e.g.  Airy’s  stress  function  in  2D.  Important 
micromechanics  observations,  that  interfacial  stress  concentrations  depend  on  the  shape  of  the  heterogeneity, 
have  been  incorporated  in  the  choice  of  stress  functions  in  (37,  88]  through  the  decomposition  of  the  stress 
functions  into  (a)  a  purely  polynomial  function  and  (b)  a  reciprocal  function  i,Trncc  (<bm  =  it™, y  +<&™c). 
The  inclusion  stress  functions  are  admitted  as  polynomial  function  <I>£0;y  ($c  =  $£oiy).  The  pure  polynomial 
function  &poiy  accommodates  the  far  field  stress  in  the  matrix  and  inclusion  and  are  written  as: 

*?oly  =  and  =  (3.3) 

P.9  P.9 

(£,  rj)  are  scaled  local  coordinates  with  origin  at  the  element  centroid  ( xc,yc ),  such  that 

£  =  (x-xc)/Le,  r)  =  (y  -yc)/Le  (3.4) 

with  the  scaling  parameter  for  each  element  Le  =  y/ max(a:  —  xc)  x  max(j/  —  yc )  V(x,  y)  G  d0.e.  The  use  of 
the  local  coordinates  (£,  if)  instead  of  global  coordinates  ( x ,  y)  in  the  construction  of  stress  functions  prevents 
ill  conditioning  of  the  stiffness  matrix  incurred  through  discrepancies  due  to  high  exponents  of  ( x ,  y)  in 
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and  <J)C.  The  reciprocal  stress  function  <l,™c  facilitates  interfacial  stress  concentration,  accounting  for  its 
shape,  and  decays  with  increasing  distance  from  the  interface. 
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The  radial  function  f(x,y)  is  a  constructed  by  Schwarz-Christoffel  conformal  transformation  of  an  elliptical 
interface  or  by  a  Fourier  series  transformation  for  arbitrary  shapes  such  as  squares  (see  [37,  68]),  to  yield 
the  properties 

f(x,  y)  =  1  on  dQ™  and  TT — V  — * ►  0  as  (2,  y)  — *  00  (3.6) 

/(z,y) 

Stress  increments  in  the  matrix  and  inclusion  phases  of  the  Voronoi  cell  elements  are  obtained  from  the 
second  derivatives  of  stress  functions  with  respect  to  x  and  y  coordinates,  resulting  in: 
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A(Txx 

A<y 

>  =  [Pm]{A/3m}  and  < 

Aaiy 

[PC]{A/3C} 


(3.7) 


where  [Pm]  and  [Pc]  are  stress  interpolation  matrices  in  the  matrix  and  inclusion.  The  boundary  displace¬ 
ments  are  generated  by  interpolation  in  terms  of  nodal  displacements  on  <9fie,  dfl™  and  d£lcc  using  linear  or 
quadratic  shape  functions. 


{Au}  =  [Le]{Aq}  on  dnc  ,  {Au"1}  =  [Lc]{Aqm}  on  (3.8) 

and  { Auc}  =  [Lc]{ Aqc}  on  dSlcc 


where  {Aq},  {Aqm}  and  {Aqc}  are  generalized  displacement  vectors  at  the  nodes.  After  substituting  stress 
and  displacement  interpolations,  the  stationary  conditions  of  equations  (3.1)  and  (3.2)  with  respect  to  the 
stress  parameters  A/3m  and  A/3C,  and  displacement  parameters  {Aq},  {Aqm}  and  {Aqc}  are  determined 
and  solved  to  yield  the  stress  and  displacement  solutions  in  each  element.  It  is  essential  that  the  relation 
between  the  interfacial  traction  components  (T^T/71)  and  separation  (um  —  uc)  be  defined  prior  to  the  above 
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step.  Different  cohesive  zone  models  are  considered  for  this  purpose. 


3.2.2  Cohesive  zone  models  for  the  interface 

Cohesive  zone  interfacial  models  [29,  6,  18]  are  effective  in  depicting  material  failure  as  a  separation  process 
across  an  extended  crack  tip.  They  introduce  softening  constitutive  equations  relating  crack  surface  tractions 
to  the  material  separation  across  the  crack.  The  tractions  across  the  interface  reach  a  maximum,  subsequently 
decrease  and  eventually  vanish  with  increasing  interfacial  separation.  A  wide  variety  of  cohesive  zone  models 
polynomial  functions  [71,  72,  73],  exponential  functions  [72,  42,  76],  bilinear  functions  [35],  and  others 
have  been  proposed  in  literature.  Motivated  by  inter-atomic  potentials  in  atomistic  modeling,  many  of 
these  cohesive  laws  use  a  potential  function  <j>  to  describe  the  traction-displacement  relation  during  material 
separation.  A  comprehensive  discussion  of  these  models  is  presented  in  [18].  Some  of  these  rate- independent 
cohesive  zone  models  are  incorporated  in  VCFEM  to  simulate  the  initiation  and  propagation  of  interfacial 
debonding. 

Polynomial  models 

Polynomial  forms  of  the  potential  functions  have  been  proposed  in  [71,  72,  73]  and  expressed  as 

if  \  27  A  Un  ,2r  4  uU\  I 

^  crmax5  {g(  )  [1  3(5*  )  +  2^  ,5*  )  1  + 

ia(|i)2[l  -  2(£)  +  (^)2]}  Vun  <  S\  (3.9) 

where  un  and  ut  are  the  normal  and  tangential  components  of  displacement  jump  across  the  interface.  In 
this  model,  the  cohesive  law  parameters  to  be  determined  from  experiments  are  omax  the  maximum  normal 
traction  in  normal  loading,  <5*,  and  a  is  the  ratio  of  interface  shear  to  normal  stiffness  at  the  interface.  The 
normal  and  tangential  components  of  the  traction  at  the  interface  are  obtained  by  taking  the  derivatives  of 
(j){un,ut)  with  respect  to  normal  and  tangential  displacement  components  as 
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Tn  = 
Tt  = 


8<f> 

dun 

d<f) 

dut 


27 


27  f  nfUn\  I  fUn\2 

-  .  crmax{a(  )[1  2(  )  +  (  ) 


(3.10) 


The  maximum  value  of  Tn  is  reached  at  8C  —  <5*/ 3.  The  traction  components  Tn  and  Tt  go  to  zero  for 
un  >  8* .  Furthermore,  the  normal  displacement  component  is  required  constrained  to  the  inequality  un  >  0 
to  avoid  penetration  between  the  fiber  and  matrix  phases  at  the  interface. 


Exponential  Model 

Cohesive  failure  models,  based  on  exponential  representation  of  cj>  in  terms  of  interfacial  separation,  have 
been  introduced  in  [76],  In  2-D,  the  model  takes  the  form 


d<j>(un,ut)  t 

J-n  —  ^  —  r  un 

UXlyX  O 

^  _  d<t>{un,ut)  _  t  n2 
h - 5 -  —  tP  ut 

OUt  o 


(3.11) 


where  8  =  \/P2u2  +  u 2  is  an  effective  displacement  and  P  assigns  different  weights  to  the  tangential  and  nor¬ 
mal  discontinuous  displacements.  The  magnitude  of  the  traction  t  =  \JT2  +  P~2T. 2  can  be  correspondingly 
expressed  as 


t  = 


(3.12) 


where  ac  is  the  maximum  cohesive  normal  traction  and  5C  is  a  characteristic  separation  length.  The  expo¬ 
nential  function  does  not  exactly  reach  zero,  and  interface  debonding  is  assumed  when  8  >  5<5C.  Unloading 
in  this  model  is  towards  the  origin,  following  a  linear  path  as 


t  =  h^-8  V<5  <  8max 

Omax 


(3.13) 
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Bilinear  model 


In  the  bilinear  model  [76],  t  is  expressed  as  a  bilinear  function  of  the  interfacial  separation  as 

t  =  <  6C  (Tmai V(5  >  6C 

dc  Oc  oe 


from  which  the  normal  and  tangential  tractions  are  derived  as 


rp  9<t> 

"  dSn  ' 

d(f)  OS  &max  r 

“  85  85n  5C  n 

VI 

Gmax  8  X 

5  5c-5en 

if  8C  <  5  <  5e 

(3.14) 

0 

if  5  >  5e 

£3 

II 

■H-©- 

ii 

d(f>  38  O max  *2  r 

85  85t  6C  P  ‘ 

if  5  <  8C 

Gmax  8  —  Se  2  r 

5  5c-s/St 

if  5C<  5  <5e 

(3.15) 

0 

if  5  >  5e 

(3.16) 

Figure  3.2(a,b)  show  the  normal  traction-separation  response  for  St  —  0  and  tangential  traction-separation 
response  for  Sn  —  0  respectively.  When  the  normal  displacement  5n  is  positive,  the  traction  at  the  interface 
increases  linearly  to  a  maximum  value  of  <jma.x  (point  A)  corresponding  to  a  value  of  6C  before  it  starts 
decreasing  to  zero  at  a  value  of  Se  (point  C).  The  unloading  behavior  in  the  hardening  region  is  linear 
following  the  loading  path.  In  the  softening  region,  the  unloading  proceeds  along  a  different  linear  path  from 
the  current  position  to  the  origin  with  a  reduced  stiffness.  This  is  shown  by  the  line  BO  in  figure  3.2(b),  for 
which  the  t  —  5  relation  is 


,  Gmax  8max  5e  c 

i  =  t - Y~Z  x  5 

Umax  Oc 


Sr  8 TTIQX  ^  dp  and  8  dn 


Reloading  follows  the  path  OBC  demonstrating  the  irreversible  nature  of  the  damage  process.  Both  normal 
and  tangential  tractions  vanish  when  8  >  Se.  The  magnitudes  of  the  tangential  traction-displacement  relation 
are  independent  of  the  sign,  and  hence  the  behavior  is  same  for  5t  positive  and  negative.  When  the  normal 
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displacement  is  negative  in  compression,  stiff  penalty  springs  with  high  stiffnesses  are  introduced  between 
the  node-pairs  at  the  interface.  Unlike  the  polynomial  and  exponential  models,  the  location  of  the  separation 
at  debonding  point  is  independent  of  the  location  of  the  peak  of  the  curve  for  the  bilinear  model.  This  gives 
flexibility  to  adjust  interfacial  parameters  for  the  peak  and  debonding  locations  to  match  the  experimental 
observations. 


3.2.3  Stress  and  Displacement  Equations  and  Solution  Methods 

Progressive  debonding  in  composite  microstructures  is  solved  using  an  incremental  approach.  In  each  in¬ 
crement,  a  set  of  element  and  global  equations  are  solved  for  stresses  and  displacements  using  the  following 
steps. 


Element  level  kinematic  equations 

Local  equations  in  each  element  are  obtained  by  substituting  stress  interpolations  (3.7)  and  displacement 
interpolations  (3.9)  in  the  element  energy  functional  of  equation  (3.1),  and  setting  variations  with  respect 
to  the  stress  coefficients  A/3m,  A/3C  respectively  to  zero.  This  results  in  the  weak  forms  of  the  element 
kinematic  relations. 


fnm  [pm]r[Sm][Pm]dfi 


[0] 


[0]  Jo  [Pc]T[Sc][Pc]dft 


/ 3m  +  A/3m 
/3C  +  A/3C 


V[pm]TKHLeI^  ~fa  nJpmiT[nc][Lc]ddfi  [°1 

[0]  [0]  /anc[Pc)T[nc][Lc]d5n 

1 

q  + Aq 


qm  +  Aqm  ) 
qc  +  Aqc 


fn,n  [pm!T{em}dfi 

Jo  [P'fVK2 


or  in  a  condensed  form 


(3.18) 


[H e]{/3  +  A/3}  =  [Ge]{q  +  Aq}  -  {R?} 


(3.19) 
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Here  [ne]  and  [nc]  are  matrices  defined  in  terms  of  direction  cosines  of  unit  outward  normal  vectors  to  the 
element  boundary  and  matrix-inclusion  interface  respectively.  The  equation  (6.42)  is  linear  and  is  solved  to 
express  the  stress  coefficients  in  terms  of  the  nodal  displacements. 


Global  traction  reciprocity  equations 


The  weak  forms  of  the  global  traction  continuity  conditions  are  subsequently  solved  by  setting  the  variation 
of  the  total  energy  functional  in  equation  (3.2)  with  respect  to  Aq,  Aqm  and  Aqc,  to  zero.  This  results  in 
the  weak  form  of  the  traction  reciprocity  conditions  as: 


N 


E 


Sanc  [Le}r[n*]T[Pm]dc>n  [0] 

~  laci^  [Lc]T[nc]T[pm]d3fl  [0] 

[o]  /anjL‘FKnp<=]dn 


P  m  +  A/3  m 

pc+Apc 


>  = 


N 


E 


/rjLe]T{*  +  At}df2 

-  fdfrm  (LC1T  ({nc}T™(u,i  +  &un,ut  +  A  ut)  +  {tc}Ttm(u„  +  Aun,  ut  +  Alii))  ddn  ’ 

-  /anc[Lc]T  ({nc}T™(un  +  A un,  ut  +  A ut)  +  {t c}T^(un  +  A un,  ut  +  A ut))  ddfl 

c  j 


(3.20) 


or  in  a  condensed  form 


N 


N 


E[Ge]T{/3  +  A/3}  =  E(R2} 


(3.21) 


Substituting  (6.42)  in  (6.44)  yields: 


N 


N 


^[GenHe]-1([G1{q  +  Aq}  -  {R?})  =  ^{Rl} 


(3.22) 
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The  normal  and  tangential  components  of  the  interfacial  separation  are  expressed  as: 


un  +  Aun  =  {nc}T[Lc]{qm  +  Aqm  —  qc  -  Aqc} 

Ut+Aut  =  {tc}T[Lc]{qm  +  Aqm-qc-  Aqc}  (3.23) 


Following  the  evaluation  of  nodal  displacements,  stress  coefficients  are  calculated  in  each  element  using  the  re¬ 
lations  (6.42).  The  stresses  at  any  location  within  the  element  may  then  be  assessed  from  the  equations  (3.7). 


Solution  Method:  The  equation  (6.45)  is  nonlinear  due  to  the  relation  between  interfacial  tractions  and 
interfacial  displacements  in  the  cohesive  laws.  A  Newton-Raphson  iteration  method  is  consequently  invoked 
to  solve  for  the  increments  of  nodal  displacement  on  the  element  boundaries  and  matrix-inclusion  interfaces. 
The  linearized  form  of  equation  (6.45)  for  the  j-th  iteration  is 


£lKT' 


dq 

dqm 

dqc 


N 


N 


e—l 


e=l 


=  -  £(GC)T«Hr1[Ge]{q  +  AqF  -  {r?}) 

=  1 

or  [Ksp{dqp  =  -  {R?ntP 


(3.24) 


This  is  iteratively  solved  to  obtain  the  incremental  nodal  displacements 


{Aqp+1  =  {Aq}>  +  {dq}>  ,  {Aqmp  +  1  =  {Aqmp  +  {dqmp  , 

{Aqcp  +  1  =  {Aqcp  +  {dqc}^  (3.25) 

The  localized  softening  in  interfacial  decohesion  can  sometimes  give  rise  to  numerical  instabilities  in  the 
Newton-Raphson  iteration,  which  is  based  on  the  smooth,  invertible  and  well-conditioned  Jacobian  iteration 
matrix,  due  to  zero  or  negative  stiffness.  The  arc-length  solver  has  been  proposed  in  [23,  24]  as  a  method 
of  overcoming  this  by  introducing  an  arc  length  as  a  replacement  to  the  incremental  load  as  the  incremental 
parameter  and  improving  the  convergence  direction  in  the  solution  space.  In  the  arc  length  method,  the 
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system  equations(3.24)  is  modified  with  the  introduction  of  the  unknown  loading  parameter  A  as 

[K 9] W}j  =  (Xj  +  dV){R9ext}  -  {RfnJJ  (3.26) 

where  dX3  and  {dq9}3  are  unknowns.  Additionally  an  orthogonality  condition  is  imposed  as 

{dq5P  •  (Aq9}J'  =  0.  (3.27) 

The  corresponding  coupled  equation  system  to  be  solved  is 

[K9]  -{RfxJ 
{Aq9P'  [0] 

Despite  advantages  of  the  arc-length  method,  the  numerical  algorithm  for  displacement  solution,  can  in 
some  cases  give  rise  to  oscillations,  especially  around  the  peak  traction  in  the  cohesive  models.  The  iterative 
solutions  do  not  converge,  with  the  slope  of  the  traction-displacement  law,  oscillating  between  large  positive 
and  negative  values.  A  regularization  method,  in  which  the  Jacobian  matrix  is  evaluated  based  on  the 
average  of  the  positive  and  negative  slopes  of  the  cohesive  models  implemented  in  this  case.  To  demonstrate 
the  effectiveness  of  this  algorithm,  an  example  with  a  square  microstructure  containing  a  single  circular 
fiber  with  a  debonding  interface  is  considered.  The  interface  uses  the  bilinear  model  in  equation(6.29)  with 
the  properties  amax  =  0.003,  5C  =  0.00002,  6e  =  0.00016,  P  =  0.707.  The  averaged  stress-strain  response 
for  the  damaging  composite  is  illustrated  in  figure  3.3.  The  Newton-Raphson  iterative  solver  stops  near 
the  peak,  but  the  arc-length  solver  continues  for  the  entire  process.  The  drop  in  stress  corresponds  to  the 
ongoing  debonding  process,  during  which  the  regular  NR  solver  is  unstable.  The  resuming  hardening  process 
corresponds  to  the  arrest  of  debonding. 


{ dq9}j 
dXj 


>  =  < 


aj{R4J  -  {RLK 
0 


(3.28) 


26 


3.2.4  Evaluation  of  volume  averaged  stresses  and  strains  for  debonding 

The  effective  macroscopic  response  of  a  composite  comprised  of  continuous  fiber  reinforcement  is  important 
to  study  the  homogenized  material  properties.  This  may  be  calculated  volume  averaging  the  local  stress  and 
strain  fields  over  the  entire  microscopic  domain  fl  (see  [61])  as 

^  J  aij(xk,t)dV 

eij(t)  =  i  /  aj(xk,  t)dV  -  (*ij(t).  (3.29) 

J  n 

where  xk  and  t  are  the  spatial  coordinates  and  time  respectively,  and 

au(0  =  f  (M01nj  +  [uj(t)]rii)ddCl  (3.30) 

^ 1  JdCl p 

oiij(t)  represents  the  effective  strain  field  caused  by  the  possible  displacement  jump  at  the  interface  due  to 
debonding.  It  is  calculated  along  the  interface  with  [u»(t)]  denoting  the  displacement  jump. 

3.2.5  Stability  of  the  VCFEM  solutions 

As  discussed  in  (37,  68,  88],  invariance  of  stresses  with  respect  to  coordinate  transformations  can  be  ensured 
by  a  complete  polynomial  representation  of  the  stress  function  <h,;.  The  necessary  conditions  for  stability 
are  that  the  tangent  compliance  modulus  in  the  strain  energy  term  be  positive  definite  and  that  the  finite¬ 
dimensional  stress  space  be  spanned  uniquely  by  the  basis  functions  ]Pm)]  and  [Pc] .  Stability  conditions  of 
the  multi-field  variational  problem  in  VCFEM  have  been  developed  in  detail  in  [88]  and  will  not  be  repeated 
in  this  chapter. 

Zhong  and  Knauss  ([117])  have  developed  a  relation  between  numerical  stability  and  material  properties 
of  the  matrix,  inclusion  and  interface  for  displacement  controlled  problems  in  ID.  For  stability,  they  have 
shown  that  the  slope  of  the  softening  region  should  be  bounded  by  a  number  that  is  determined  in  terms 
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material  and  interface  properties  and  the  size  of  the  body  L,  or  mathematically 


2  E  dT(S) 

^  T71CLX[  —  ^  ]| softening 


(3.31) 


where  T{5)  is  the  cohesive  traction  and  E  =  Ei  and  E2  being  the  Young’s  modulus  of  the  two 

materials  bonded  by  the  cohesive  zone.  Larger  softening  slopes  in  the  cohesive  model  make  this  criterion 
difficult  to  satisfy.  To  examine  this  characteristic,  a  1  x  1  square  domain  with  a  circular  fiber  of  radius 
0.5,  is  simulated  for  tension  loading.  The  interface  is  lined  with  the  polynomial,  exponential  and  two 
bilinear  cohesive  models  as  discussed  in  section  3.2.2.  The  material  properties  for  matrix  and  fiber  are 
Ematrix  =  4.6  GPa,  v matrix  =  0.4,  Einciusion  =  210  GPa,  ^inclusion  =  0.3  respectively.  The  cohesive  zone 
model  parameters  for  all  models  are  chosen  to  satisfy  the  stability  condition  stabcrit,  and  are  tabulated 
in  table  3.5.  For  the  same  ac  and  6C,  the  polynomial  model  has  a  faster  rate  of  decay  to  zero  («  2<5C) 
in  comparison  with  the  exponential  model  (?»  4 8C).  The  bilinear  models  are  chosen  with  Se  =  3 5C  in  (1) 
and  Se  =  5<5C  in  (2).  The  corresponding  macroscopic  stress-strain  responses  are  illustrated  in  figure  3  4. 
All  models  are  able  to  simulate  the  entire  debonding  process  from  hardening  to  post-debonding.  However, 
the  polynomial  model  with  the  largest  softening  slope  undergoes  a  sudden  drop  near  the  peak  and  some 
oscillations,  before  it  stabilizes  again. 

3.2.6  Adaptive  Enhancement  of  the  Voronoi  Cell  FE  Model 

To  establish  rapid  convergence  of  microscale  solutions  and  enhance  the  solution  accuracy,  a-posteriori  adap¬ 
tivity  has  been  incorporated  in  the  Voronoi  cell  finite  element  model  without  debonding  in  [88].  Adaptations 
based  on  suitably  chosen  error  indicators  are  introduced  as  follows. 

i.  Adaptation  to  reduce  traction  reciprocity  error  on  element  boundaries  and  interfaces:  To  estimate  the  quality 
of  solution  induced  by  the  weak  satisfaction  of  traction  continuity  on  Voronoi  cell  element  boundary,  an 


28 


average  traction  continuity  error  ( A.T.R.E .)  is  defined  as: 


A.T.R.E. 


2^e=l 


„e  , 

-T  +  2~,c=\ 


Ne  +  Nc 


(3.32) 


where 


e 


C 

r 


1  Jan.CtW]  • 

a  n *  fdQ  ddfle 


and 


°  n*  fan,  mi 


(3.33) 


In  equation  (3.32)  Ne  and  Nc  are  the  total  of  all  segments  on  all  element  boundaries  dQ.e  and  interfaces 
SHc  respectively.  The  stress  d  in  the  denominator  is  the  absolute  maximum  principal  value  of  the  volume 
averaged  stress  tensor  in  the  microstructure,  viz.  dij  =  and  n*  the  number  of  degrees  of  freedom  per 

node  in  the  problem  and  [|t|j  is  the  traction  discontinuity  along  different  element  boundaries  and  interfaces 
in  the  microstructural  model.  The  traction  continuity  error  in  equation  (3.32)  is  minimized  by  selectively 
enhancing  boundary  and  interface  displacement  degrees  of  freedom  in  the  directions  of  optimal  displacement 
enrichments.  These  directions  minimize  the  virtual  work  due  to  traction  discontinuity  and  are  obtained  from 
components  of  the  traction  discontinuity  in  directions  orthogonal  to  the  original  displacement  field. 


ii.  Adaptation  to  improve  stress  concentration  at  the  crack  tip  on  the  interface:  As  seen  in  figure  3.5(a), 
node  pairs  are  initially  positioned  at  equal  arc-lengths  along  the  interface.  Adaptation  for  reducing  error  in 
traction  continuity  puts  additional  nodes  on  the  interface  as  discussed  in  item  (i).  The  node  positioning  is 
thus  far  independent  of  the  extent  of  debonding.  However,  the  crack  tip  stresses  are  better  represented  if  the 
nodes  are  coincident  or  at  least  near  the  crack  tip.  Correspondingly,  a  set  of  nodes  in  this  model  are  moved 
with  the  evolution  of  the  cohesive  zone  to  provide  high  resolution  in  the  regions  of  high  cohesive  tractions 
across  interfaces  during  debonding.  As  seen  in  figure  3.5(a),  when  the  peak  of  the  cohesive  zone  model  lies 
between  two  neighboring  nodes,  one  or  both  of  them  are  moved  close  to  the  peak  traction.  For  two  adjacent 
node  pairs  n  and  n  +  1,  if: 

5n  >  5C  and  (5n.fi  <  Sc  or  (5n  <  5C  and.  <5n+i  >  dc  (3.34) 
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then  the  critical  point  with  a  displacement  jump  5C  exists  between  the  two  nodes.  For  a  linear  mapping,  the 
coordinates  of  the  critical  point  are 


Xc  —  Xn  "1" 


5C  -  Sn 

5  71  f  1  5  7 


yc  =  yn  + 


5C  -  5n 

^n+1  5  ji 


(l/n+1  2/n) 


(3.35) 


The  node  n  may  be  moved  to  this  location  to  generate  an  optimal  stress  representation.  A  similar  interpo¬ 
lation  can  be  applied  to  interpolate  the  displacements  at  the  critical  point  from  the  n  —  th  and  n  +  l  —  th 
nodes.  The  results  of  the  adaptation  by  the  above  two  techniques  are  illustrated  in  figure  3.5(b),  for  the 
square  matrix  with  a  single  circular  fiber.  The  maximum  error  in  the  traction  reciprocity  as  a  function  of  the 
macroscopic  strain  is  plotted  in  these  figures.  The  error  on  the  element  boundary  is  reduced  considerably 
by  node  adaptivity  that  increases  the  number  of  nodes  from  four  to  fifteen  on  the  element  boundary.  The 
other  set  of  results  correspond  to  reduction  in  traction  error  on  the  interface  due  to  node  movement  along 
the  interface.  The  effectiveness  of  the  adaptation  techniques  is  adequately  established  in  these  examples. 


3.3  Numerical  Examples 

3.3.1  Model  Validation 

Prior  to  its  application  in  problems  of  multi-fiber  composite  microstructures,  the  VCFE  model  with  the 
cohesive  zone  interface  is  validated  for  accuracy  and  reliability.  In  the  first  example  VCFEM  results  are 
compared  with  those  in  [69]  for  interfacial  crack  initiation  and  growth  in  a  transversely  loaded  composite. 
The  microstructure  is  represented  by  a  unit  cell  in  a  uniform  hexagonal  array  with  a  fiber  volume  fraction 
Vf  =  0.5,  as  shown  in  figure  3.6(a).  The  unit  cell  dimensions  are  shown  in  the  figure  where  r  =  10 gm  and 
b  =  15.55 fim.  Isotropic  and  linear  elastic  matrix  and  fiber  properties  are  given  in  table  3.5.  The  matrix-fiber 
interface  is  represented  by  linear  elastic  springs  prior  to  failure.  The  normal  and  tangential  tractions  are 
assumed  to  be  independent  and  are  proportional  to  displacement  jumps  in  the  two  directions  respectively, 
i.e. : 

Tn  =  M«n]  .  Tt  -  kt[ut]  V  [w«l  ^  0  (3.36) 
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where  [.]  refers  to  jump  across  the  interface,  and  kn  and  kt  are  normal  and  tangential  stiffness  constants. 
The  impenetrability  constraint  [u„]  >  0  is  enforced  using  a  penalty  spring.  Debonding  is  assumed  to  initiate 
in  this  model  when  the  normal  displacement  jump  [un]  reaches  a  critical  value  5  =  lOnm.  Post-debonding 
tractions  are  assumed  to  be  zero,  i.e. 


Tn  =  Tt  =  0  for  un>  6  (3.37) 

The  macroscopic  stress-strain  response,  generated  by  VCFEM  is  compared  with  that  in  [69]  in  figure  3.6(b). 
The  two  curves  agree  very  well  with  a  small  difference  during  debonding. 

In  the  second  example,  VCFEM  predictions  are  compared  with  experiments  on  debonding  of  composites 
that  have  been  described  in  [37].  The  experiment  simulated  here  is  conducted  with  a  single  specimens 
in  the  form  of  a  cruciform  as  shown  in  figure  3.7(a).  The  cruciform  shape  has  been  developed  to  avoid 
stress  singularity  at  the  intersection  of  fiber-matrix  interface  and  free  surface,  that  occur  in  uniform  width 
specimens.  The  most  significant  advantage  of  the  cruciform  geometry  is  that  it  forces  debond  failure  to 
initiate  in  the  central  region  of  the  specimen.  The  reinforcing  fibers  in  the  composite  system  are  stainless 
steel  filaments,  while  the  matrix  is  an  epoxy  resin  that  is  cured  with  polyetheramine.  The  epoxy  matrix  is 
transparent  and  allows  visualization  of  the  debonding  process  at  the  fiber-matrix  interface.  A  very  thin  film 
of  weak  strength  freekote  (<  0.1  fim)  is  inserted  as  the  interface  material.  This  allows  a  somewhat  stable 
growth  of  the  debond  crack.  The  model  specimens  are  loaded  in  tension  on  a  servo-hydraulic  testing  machine. 
The  onset  of  fiber-matrix  debonding  is  identified  with  a  sharp  change  in  the  slope  of  the  experimental  stress- 
strain  curve  in  figure  3.7(b).  Subsequent  loading  proceeds  with  a  lower  stress-strain  slope,  due  to  a  reduced 
load  carrying  capability  of  the  partially  debonded  fiber.  Following  failure,  the  specimens  are  sectioned  at 
the  center  and  fluorescent  dye  penetration  is  used  to  determine  the  total  angle  of  debond  as  approximately 
85°  as  shown  in  figure  3.7(c).  The  cohesive  zone  parameters  have  been  determined  in  [37]  by  solving  an 
inverse  problem  in  which  the  difference  between  the  experimental  and  simulated  results  is  minimized.  For 
the  bilinear  model  the  cohesive  parameters  are:  <7max  =  0.0037 GPa,  5C  =  0.0028,  5e  =  0.0035,  and  fi  —  0.707. 
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The  geometric  properties  of  the  composite  simulated  include  a  6.82  mm  x  6.0  mm  domain  with  a  circular 
fiber  of  radius  2.36  mm.  The  material  properties  are:  Esteei  =  210 GPa,  vsteei  =  0.3,  Eepoxy  =  4.6 GPa  and 
Vepoxy  =  0.4.  The  macroscopic  stress-strain  plot  of  the  simulation  in  figure  3.7(b),  match  the  experimental 
result  very  well.  Furthermore,  the  simulated  debonding  angle  in  the  composite  microstructure  is  obtained 
to  be  90°  as  shown  in  figure  3.7(d).  Both  macroscopic  and  microscopic  results  obtained  from  VCFEM 
simulations  are  found  to  yield  satisfactory  comparison  with  experiments. 

3.3.2  Effect  of  interfacial  properties  on  debonding 

The  shape  of  the  traction  displacement  curve  in  the  cohesive  zone  model  plays  an  important  role  in  the 
simulation  of  initiation  and  progress  of  debonding  in  composite  microstructures.  The  total  cohesive  energy 
has  been  split  into  an  intrinsic  energy  dissipation  (rtni)  and  an  extrinsic  cohesive  energy  dissipation(reIt) 
in  [18],  depending  on  the  ascending  and  descending  portions  of  the  curve.  The  ratio  of  the  extrinsic  and 
intrinsic  cohesive  energy  energies  ip  =  is  used  to  denote  a  shape  factor  for  the  cohesive  law.  The 
VCFEM  with  the  bilinear  cohesive  law  is  used  to  study  debonding  in  (i)  a  single  fiber  microstructure,  (ii) 
a  random  microstructure  with  a  cluster,  and  (iii)  a  real  micrograph  of  fiber  reinforced  composite  with  264 
fibers. 

Microstructure  with  a  single  fiber 

The  single  fiber  microstructure  of  the  previous  section  with  the  bilinear  cohesive  model  interface,  is  simulated 
with  four  different  sets  of  cohesive  zone  parameters  listed  in  table  3.5.  As  shown  in  the  inset  of  figure  3.8, 
the  sets  A,  B  and  C  have  the  same  peak  stress  (<xmaI),  whereas  the  set  D  has  a  peak  stress  of  2amax.  The 
shape  factor  ip  for  sets  B  and  D  are  2.0,  while  those  for  A  and  C  are  smaller.  The  macroscopic  response 
for  the  different  interface  laws  are  illustrated  in  figure  3.8.  The  post-debonding  responses  are  almost  same 
for  all  cases.  However,  significant  differences  exist  in  the  softening  region.  The  case  B  with  larger  ip  and 
smaller  stiffness  (due  to  the  softer  interface),  compared  to  cases  A  and  C debonds  later.  The  case  D  with  a 
larger  amax  exhibits  a  totally  different  macroscopic  response  in  the  softening  region.  The  initial  slope  of  a 
macroscopic  stress-strain  plot  depends  on  the  hardening  slope  of  the  cohesive  law,  because  of  stress  transfer 
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at  the  interface.  The  subsequent  failure  behavior  is  determined  by  the  softening  part  of  the  cohesive  curve. 
For  the  same  peak  amax  and  the  total  cohesive  energy,  the  debonding  behavior  has  a  considerable  dependence 
of  the  shape  factor  ip.  Larger  drops  in  the  stress-strain  plots  are  observed  for  larger  shape  factors. 

Random  distributions  with  and  without  cluster 

This  example  is  constructed  to  examine  effect  of  interfacial  laws  on  the  debonding  behavior  of  multi-fiber 
microstructures.  Specifically,  a  random  microstructure  and  a  clustered  microstructure,  both  with  49  fibers 
as  shown  in  figures  3.9(a)  and  (b)  are  considered. 

The  cluster  contains  8  fibers.  Three  sets  of  cohesive  law  parameters  are  implemented  for  each  microstruc¬ 
ture  as  listed  in  table  3.5.  The  sets  A  and  C  have  the  same  crmax  and  6e,  but  the  shape  factor  ip  for  C 
is  larger.  The  peak  stress  <7max  for  B  is  relatively  smaller.  As  for  the  single  fiber  case,  the  macroscopic 
stress-strain  plots  in  figures  3.9(c)  and  (d)  show  considerable  dependence  on  the  shape  factor  and  peak 
stress.  Damage  initiation  and  propagation  in  the  microstructure  is  affected  considerably  by  the  shape  of  the 
cohesive  zone  models  for  two  spatial  distributions.  The  contour  plots  of  axx  in  figure  3.10  illustrate  different 
debonding  propagation  paths  for  the  different  interfaces  and  distributions. 

The  locations  of  the  debond  initiation  is  the  same  for  the  different  cohesive  laws.  For  the  clustered 
distribution,  the  case  C  shows  a  clear  debond  propagation  path  in  comparison  with  the  other  two. 

3.3.3  A  real  micrograph  with  264  fibers 

In  this  example,  an  optical  micrograph  of  a  real  composite  shown  in  figure  3.11(a),  is  modeled  by  VCFEM. 
The  micrograph  is  adapted  to  a  computational  domain  and  tessellated  into  a  mesh  of  Voronoi  elements, 
as  shown  in  figure  3.11(b).  The  material  properties  for  matrix,  fiber  and  interface  are:  Esteei  —  210 GPa, 
lasted  =  0.3,  Ecpoxy  =  4.6GPa,  vepoXy  =  0.4,  ac  =  5 MPa,  <5C  =0.000051,  and  (3— 0.707.  Simulations  with  two 
different  cohesive  law  shape  factors,  ip\  =  0.138  and  ip2  =  0.197,  generate  two  different  debonding  paths.  A 
well  defined  damage  path  is  observed  in  figure  3.11(d)  for  the  case  with  the  higher  shape  factor.  For  the 
lower  ip i ,  a  bifurcation  is  seen  in  the  damage  paths. 
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The  above  examples  show  conclusively  that  in  addition  to  the  cohesive  energy,  amax  and  6e,  the  shape 
factor  plays  an  important  role  in  determining  the  macroscopic  softening  response  as  well  as  the  debonding 
path.  When  ip  is  large  the  debonding  occurs  later  and  there  is  a  sudden  drop  in  the  macroscopic  stress-strain 
response.  Furthermore,  the  path  of  microscopic  debond  propagation  is  more  defined  for  higher  shape  factors. 

3.4  Microstructural  Characteristics  on  the  Initiation  of  Debond¬ 
ing 

From  a  microstructure  design  perspective,  two  aspects  are  of  interest  to  the  composite  community.  The  first 
aspect  is,  at  what  macroscopic  strain  does  microstructural  debonding  initiate  for  a  given  multi-fiber  config¬ 
uration.  The  second  is  the  location  where  this  occurs.  In  other  words,  what  are  the  local  microstructural 
characteristics  that  trigger  interface  failure.  Microstructural  characterization  of  non-homogeneous  compos¬ 
ites  has  been  conducted  in  [39,  40]  using  various  statistical  functions  of  geometric  parameters.  These  include 
the  cumulative  function  and  probability  density  functions  of  local  area  fraction  and  near-neighbor  distances, 
the  second  order  intensity  function  K(r)  and  the  pair  distribution  function  g(r)  etc.  These  functions  help 
identify  spatial  distributions  like  uniform,  random  or  clustered  patterns.  The  present  example  is  aimed  at 
the  study  of  the  effect  of  microstructural  morphology  on  the  damage  evolution. 

A  special  function  is  constructed  in  this  example  as  the  weighted  sum  of  various  geometric  parameters 
that  can  contribute  to  the  initiation  of  debonding.  For  the  k  —  th  fiber,  the  function  is  defined  as: 

n 

Gk  =  J2^iSi  (3-38) 

«=i 

where  Wi  is  the  weight,  and  Si  are  specific  geometric  parameters  describing  the  local  spatial  distribution.  In 
this  work,  four  different  parameters  are  used  and  hence  n  =  4.  For  a  domain  including  N  fibers  or  inclusions, 
normalized  parameters  are  defined  as 
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1.  Sk\  is  a  measure  of  the  normalized  local  area  fraction  for  the  k- th  fiber. 


(. LAF)k  -  Min  ( LAF)j 

qk  _  _ *<j<N _ 

1  Max  (LAF)j  -  Min  (LAF)i 

l<j<N  l<j<N 


(3.39) 


where  N  is  the  total  number  of  fibers  and  ( LAF)i  is  the  local  area  fraction  for  the  j-th  fiber,  which 
is  evaluated  as  the  ratio  of  the  fiber  cross-sectional  area  to  the  area  of  the  respective  Voronoi  cell  (see 
[39,  40]). 

2.  Sk:  is  a  measure  of  the  normalized  inverse  of  the  nearest  neighbor  distance  for  the  &-th  fiber 


(. INND)k  -  Min  (INND)i 

nk  _  _ i<j<n _ 

2  Max  (INND)j  -  Min  (INND)i 

l<j<N  1  <j<N 


(3.40) 


where  ( INND )3  is  the  inverse  of  nearest  neighbor  distance  of  j-th  fiber.  The  near  neighbors  of  a  given 
fiber  are  those  that  share  common  edges  of  the  Voronoi  cell. 


3.  Sf:  is  a  measure  of  the  normalized  fcth  fiber  size. 

(FS)k  -  Min  (FSy 

sk  _  i<i<N _ 

3  Max  (FSV  -  Min  (FS)j 


(3.41) 


where  (FS)-7  is  the  area  of  ^-th  fiber. 

4.  S4 :  is  a  measure  of  the  normalized  average  size  of  fibers  around  the  fc-th  fiber. 

(AFS)k  -  Min  (AFSy 

qk  _  _ '<i<N _ _ 

4  Max  (AFS)i  -  Min  (AFSy 

l<j<N  1  <j<N 


(3.42) 


where  (AFS)3  is  the  average  area  of  fibers  around  j-th  fiber.  The  range  of  each  Sk  is  [0,1]  in  the  above 
definitions.  Both  S3  and  Sk  are  zero  for  microstructures  containing  same  sized  fibers,  but  they  affect 
microstructures  with  size  variations.  The  weights  in  the  equation  3.38  are  assigned  selectively  for  high  effec¬ 
tiveness  of  the  characterization  function  Gk.  Four  representative  microstructures  under  simple  tension  are 
considered  for  evaluating  the  effect  of  morphology  on  debonding.  They  are  (a)  random  microstructure  with 
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49  equi-sized  fibers  of  figure  3.9(a),  (b)  the  clustered  microstructure  with  49  equi-sized  fibers  of  figure  3.9(b), 

(c)  a  random  microstructure  with  38  variable-sized  fibers  (ratio  of  maximum  to  minimum  radius  =1.62),  and 

(d)  a  clustered  microstructure  with  55  variable-sized  fibers  (ratio  of  maximum  to  minimum  radius  =2.64). 
The  bilinear  cohesive  law  with  ac  =  12 MPa,  Sc  =  0.000052,  Se  =  0.000094,  and  fi=  0.707  is  used  for  the 
interfaces.  Material  properties  are  Efn,er  =  210 GPa,  v fiber  =  0.3,  Ematrix  —  4.6 GPa  and  vmatrix  =  0.4. 

To  understand  the  effect  of  each  of  the  parameters  Sk  on  the  initiation  of  debonding  in  the  different 
microstructures,  a  sensitivity  study  is  conducted.  In  this  study,  debonding  simulations  are  conducted  for 
each  microstructure  and  the  location  of  initiation  is  noted.  Subsequently,  the  characterization  function  Gk 
with  the  weights  wl  set  to  1  and  0  to  manifest  each  Sk.  A  summary  of  the  results  for  each  of  the  Sk's  is 
provided  in  table  3.5.  For  the  microstructures  with  equi-sized  fibers,  S 1  and  S2  are  good  indicators  of  the 
initial  debonding  location,  while  S 3  and  S 4  do  not  since  they  reflect  size  effects.  On  the  other  hand,  for 
the  microstructures  with  variable  sizes,  S3  and  5 4  are  better  predictors  of  the  initial  debonding  location  in 
comparison  with  S 1  and  S2.  Interfacial  debonding  is  thus  sensitive  to  the  fiber  size  in  addition  to  local  area 
fraction  and  nearest  neighbor  distance.  The  weights  Wi  are  adjusted  to  optimal  values  after  several  iterations 
and  are  uq=0.1,  102=0.4, 103=1.5,  and  104=1. 5.  Figure  3.5  shows  the  contour  plots  for  Gk  and  the  locations 
of  initial  debonding  with  the  stresses  in  the  load  direction.  The  geometric  indicator  is  found  to  catch  the 
initiation  location  adequately.  Finally  the  macroscopic  debond  initiation  strain  is  presented  against  values  of 
the  maximum  Gk  in  each  microstructure  in  table  3.5.  It  is  clear  that  the  function  Gk  has  a  direct  bearing  on 
the  strain,  and  larger  values  signal  a  smaller  failure  strain.  This  implies  the  strength  of  this  characterization 
function  is  predicting  microstructural  failure. 

3.5  Conclusions 

Numerical  simulations  are  conducted  with  the  polynomial,  exponential  and  bilinear  cohesive  laws  to  un¬ 
derstand  the  effect  of  cohesive  laws  on  the  debonding  process.  It  is  observed  that  in  addition  to  the  total 
cohesive  energy,  the  shape  of  the  traction-displacement  has  an  effect  on  initiation  and  especially  on  propaga- 
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tion  of  debonding  through  the  microstructure.  The  bilinear  model  which  has  more  flexibility  with  respect  to 
ascending  and  descending  portions  of  the  cohesive  curve  is  chosen  for  subsequent  debonding  simulations.  The 
sensitivity  of  debonding  in  random  and  clustered  microstructures  to  various  shapes  in  the  bilinear  cohesive 
model  is  studied.  The  clustered  microstructures  initiate  debonding  at  smaller  strains.  Also  for  the  clustered 
microstructure,  the  debonding  is  observed  to  proceed  along  defined  paths  rather  than  in  a  random  manner 
especially  for  cohesive  curve  with  larger  softening  energy.  Finally,  a  quantitative  characterization  function  is 
developed  in  terms  of  weighted  microstructural  geometric  features  like  local  area  fraction,  nearest  neighbor 
distance,  fiber  size  and  neighboring  fiber  size,  to  predict  the  location  of  damage  initiation.  The  sensitivity 
of  the  damage  to  the  various  geometric  features  is  utilized  to  determine  the  weights.  Result  show  that  the 
function  is  quite  successful  in  predicting  the  location  of  damage  onset  with  good  accuracy.  This  study  reveals 
the  significance  of  analyzing  large  regions  of  the  microstructure  and  proves  the  effectiveness  of  the  VCFEM 
analysis  for  the  same.  The  Voronoi  cells  also  provide  the  essential  link  between  the  microstructural  features 
and  response  that  is  important  in  damage  analysis. 
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Polynomial  Exponential 


Bilinear 

Parameter  I  Parameter  II 


@max  {GPcl) 

0.003 

0.003 

0.003 

0.003 

Sc 

0.000167 

0.000167 

0.000167 

0.000167 

Se 

0.0005 

0.000835 

0.0005 

0.000835 

a 

p 

10.0 

0.707 

0.707 

0.707 

Table  3.1:  Interfacial  properties  for  various  cohesive  zone  models. 


EA(GPa)  uA  ET(GPa )  uT  GA(GPa)  GT{GPa)  KT(GPa) 


Graphite  fiber  232 
Epoxy  matrix  5.35 


15.0 

0.49 

24.0 

5.35 

0.354 

1.976 

1.9 

Table  3.2:  Material  elastic  properties  for  simulation  of  hexagonal  unit  cell  in  Moran(1991). 
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A 

B 

C 

D 

&max  (GPa) 

0.003 

0.003 

0.003 

0.006 

<5C 

0.0001 

0.0002 

0.00005 

0.0001 

Sc 

0.0003 

0.0003 

0.0003 

0.00015 

P 

0.707 

0.707 

0.707 

0.707 

■0 

0.5 

2.0 

0.2 

2.0 

Table  3.3:  Cohesive  law  parameters  with  variation  in  the  shape. 


A 

Random 

B 

c 

A 

Cluster 

B 

C 

&max  (GPa) 

0.0012 

0.0006 

0.0012 

0.0012 

0.0006 

0.0012 

<5c 

0.000052 

0.000104 

0.000104 

0.000052 

0.000104 

0.000104 

<5e 

0.000145 

0.000145 

0.000145 

0.000195 

0.000195 

0.000195 

P 

0.707 

0.707 

0.707 

0.707 

0.707 

0.707 

4> 

0.559 

2.54 

2.54 

0.364 

1.18 

1.18 

Table  3.4:  Interfacial  Properties. 


Weights 

Random 

Cluster 

Variable  size 

Variable  size  with  Cluster 

U)\ 

w2 

' w3 

W4 

1 

0 

0 

0 

Strong 

Strong 

Weak 

Weak 

0 

1 

0 

0 

Strong 

Strong 

Weak 

Weak 

0 

0 

1 

0 

Weak 

Weak 

Strong 

Strong 

0 

0 

0 

1 

Weak 

Weak 

Strong 

Strong 

Table  3.5:  Sensitivity  of  debonding  initiation  to  parameter  Sf. 


Random 

Cluster 

Varying  size 

Varying  size  &  Cluster 

Maximum  Gk 

0.436 

1.21 

Strain 

0.0018 

0.0015 

Table  3.6:  Debonding  initiation  and  characterization. 
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Figure  3.3:  Macroscopic  stress-strain  response  demonstrating  the  improvement  with  arc-length  stability. 


Figure  3.4:  Comparison  of  macroscopic  stress-strain  response  for  various  cohesive  models. 
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Maximum  discontinuous  traction(MPa) 


Figure  3.5:  (a)  Node  adaptation  on  element  boundary  and  interface  (b)  maximum  traction  discontinuity  on 
element  boundary  and  interface  before  and  after  adaptation 


fa) 


(b) 


Figure  3.6:  (a)  Unit  cell  model  of  the  microstructure  with  hexagonal  array,  (b)  comparison  of  macroscopic 
stress-strain  response. 
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(a) 


(b) 
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5.248E-03 


-2.440E-03 


(c)  (d) 

Figure  3.7:  Comparison  of  simulation  with  experiment;  (a)  faceview  of  the  debonded  cruciform  specimen 
showing  dye  penetration,  (b)  comparison  of  macroscopic  stress-strain  response,  (c)  the  cross  section  indi¬ 
cating  debonding  angle  as  the  limits  of  the  dye  penetrated  region,  (d)  contour  plot  of  the  microscopic  axial 
stress(GPa) 
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Macroscopic  Stress  (GPa) 


Figure  3.8:  Macroscopic  stress-strain  response  for  various  shape  factors  of  the  cohesive  zone  models. 


3.W3E-02 


(b) 


(e) 


(0 


Figure  3.10:  Contour  plots  showing  the  stress  <7n(GPa)  and  damage  evolution  in  the  microstructure:  (a) 
random  microstructure  with  cohesive  property  A  at  e=0. 00180,  (b)  clustered  microstructure  with  cohesive 
property  A  at  e=0. 00270.  (c)  random  microstructure  with  cohesive  property  B  at  e=0. 00180.  (d)  clustered 
microstructure  with  cohesive  property  A  at  e=0. 00270.  (e)  random  microstructure  with  cohesive  property 
C  at  e=0.00195.  (f)  clustered  microstructure  with  cohesive  property  A  at  e=0. 00270. 


48 


Figure  3.11:  (a)  Optical  micrograph  of  a  real  composite  with  264  fibers,  (b)  computational  model  incorpo¬ 
rating  the  Voronoi  cell  mesh,  (c)  oxx{GV&)  contour  plot  at  6^=0.00138  for  «5e=0.00042  and  0=0.138,  (d) 
On(GPa)  contour  plot  at  e„=  0.000967  for  <5e=0. 00031, and  0=0.197. 
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Chapter  4 

The  Extended  Voronoi  Cell  Finite 
Element  Model  For  Multiple  Cohesive 
Cracks  Propagation 

4.1  Introduction 

Numerical  analysis  and  simulation  of  the  growth  of  multiple  cracks  in  materials  is  a  challenging  enterprise 
due  to  morphological  and  constitutive  complexities  that  govern  its  growth.  The  conventional  finite  element 
method  suffers  from  very  slow  convergence  since  the  element  formulation  does  not  account  for  high  gradients 
and  singularities.  Even  a  very  high  density  mesh  cannot  overcome  pathological  mesh  dependence  near  the 
crack  tips  and  avoid  biasing  the  direction  of  crack  propagation. 

In  this  chapter,  an  extended  VCFEM  or  X-VCFEM  is  developed  for  modeling  the  growth  of  multiple  cohe¬ 
sive  cracks  in  a  brittle  material.  The  model  accounts  for  interaction  between  cracks  and  invokes  an  adaptive 
crack  growth  formulation  to  represent  the  continuously  changing  direction  of  evolving  cracks.  X-VCFEM 
augments  the  conventional  VCFEM  model  by  incorporating  multi-resolution  wavelet  functions  [41,  51,  80] 
in  the  vicinity  of  the  crack  tip,  in  addition  to  branch  functions  based  on  level  set  methods.  The  incremental 
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crack  propagation  direction  and  length  are  adaptively  determined  by  a  cohesive  energy  based  criterion.  No 
remeshing  is  needed  in  X-VCFEM  for  simulating  crack  growth,  and  this  adds  to  its  desirability  and  effec¬ 
tiveness.  It  begins  with  the  X-VCFEM  formulation,  followed  by  numerical  example  showing  the  convergence 
of  this  model.  Then,  X-VCFEM  is  used  to  understand  the  influence  of  cohesive  parameters,  e.g.  peak 
stress  and  critical  separation  on  crack  growth  in  a  monolithic  brittle  material.  Subsequently,  the  effect  of 
morphological  distributions  including  crack  interaction,  clustering,  alignment,  etc.  on  growth  and  merging 
are  studied  as  important  factors  critical  to  the  failure  process. 

4.2  Voronoi  Cell  Fem  Formulation  for  Multiple  Propagating  Cracks 

The  Voronoi  cell  finite  element  mesh  for  a  brittle  matrix  with  a  dispersion  of  pre-existing  cracks  is  shown 
in  figure  4.1(a).  The  typical  Voronoi  cell  mesh  corresponds  to  an  unstructured  mesh  that  is  generated  by 
Dirichlet  or  Voronoi  tessellation  of  the  domain,  based  on  the  position,  shape  and  size  of  heterogeneities  (in¬ 
clusion,  void,  crack  etc.).  Various  tessellation  schemes  have  been  discussed  and  developed  in  [37,  68].  While 
the  name  Voronoi  cell  has  been  historically  used  because  of  its  association  with  point  seeds  in  the  generation 
process,  the  cells  used  in  VCFEM  may  be  variants  of  this  construct.  Essentially  they  represent  neighbor¬ 
hood  or  regions  of  influence  for  each  heterogeneity.  Subsequently  the  Voronoi  cell  FE  formulation  considers 
each  cell  as  a  super-element  consisting  of  a  heterogeneity  and  its  neighborhood  surrounding  matrix  [68,  88] 
without  any  further  subdivision.  The  interfacial  debonding  analyses  in  [37,  58]  invoke  the  cohesive  zone 
models  to  represent  the  growth  of  interfacial  crack.  However  the  main  difference  between  that  formulation 
and  the  present  one  is  that,  in  the  present  case  the  path  of  the  crack  is  arbitrary  and  is  a-priori  unknown. 
This  poses  significant  challenges  that  have  been  overcome  with  the  X-VCFEM  formulation. 

Consider  a  pre-cracked  microstructural  region  fl  consisting  of  N  cracks  as  shown  in  figure  4.1(a).  The 
region  is  divided  into  an  unstructured  finite  element  mesh  of  arbitrary  Voronoi  cells.  A  typical  VC  element 
Cle  containing  a  crack  and  its  neighboring  matrix  is  depicted  in  figure  4.1(b).  The  element  boundary  dfle  with 
outward  normal  nE  may  consist  of  regions  of  prescribed  traction  Tte,  prescribed  displacement  rue  and  inter- 
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element  edges  rme,  i.e.  dfl f  =  r(e(jrue  Urme.  Furthermore,  each  element  consists  of  a  crack  containing 
a  fracture  process  zone  that  is  represented  by  a  cohesive  zone  model.  The  incompatible  displacement  field 
across  the  crack  rcr  is  facilitated  through  a  set  of  connected  node- pairs  along  the  crack  length.  The  node-pair 
merges  at  the  crack  tip  by  enforcing  the  same  displacement.  The  normal  along  the  crack  path  is  denoted  by 
ncr.  For  the  VCFEM  element  formulation,  the  micromechanics  boundary  value  problem  is  described  as: 


Find  (cr,u'E,ucr)  €  T  x  VE  x  Vcr  satisfying 
d  13 

V  •  a  +  f  =  0  and  =  e  G  fle 

da 

u£  =  u  on  rue  ,  a  nE  =  t  on  F(e  and  a  ■  n”-  =  t‘:D,l  on  Tcr 


(a) 

(b) 


(4.1) 


The  variables  a,  e,  B  and  f  are  the  equilibrated  stress  fields,  the  corresponding  strain  fields,  the  complimen¬ 
tary  energy  and  body  forces  per  unit  volume  respectively  in  the  element  interior.  T,  V£  and  Vcr  correspond 
to  Hilbert  spaces  containing  the  stress  and  displacement  solutions  respectively.  u£  is  the  kinematically 
admissible  displacement  field  on  the  element  boundary  dO.E  and  ucr  represents  the  displacements  on  the 
internal  cohesive-crack  surfaces  rcr.  Variables  with  superscript  E  are  on  the  element  boundary  while  those 
with  superscripts  cr  correspond  to  the  crack  surface.  The  traction  tcoh  between  node-pairs  on  the  crack 
surface  are  modeled  by  the  cohesive  zone  traction-separation  law.  The  VCFEM  formulation  is  based  on 
the  assumed  stress  hybrid  finite  element  method,  in  which  stationarity  conditions  of  the  element  energy 
functional  in  the  variational  principle  yields  weak  forms  of  the  kinematic  equation  and  traction  reciprocity 
conditions,  as  Euler  equations.  In  the  small  deformation  elasticity  incremental  formulation  for  evolving 
cracks,  the  element  energy  functional  ne  is  defined  in  terms  of  increments  of  stresses  and  displacements  as: 


+ 

+ 


ne(ffy,  A<ry,uf,Auf,«f  .Auf)  =  -  /  AB(aij,Aaij)dQ.  -  f 

J  fic 

I  (aij  +  Aoij)nE(uE  +  Auf)ddQ.  -  f  +  Ati)(uEi  +  AuEi)drte 
Jan,  Jrtrn 

/,  (aij  +  Aaij)nf{uf  +  Auf)drcr  -  f2  ( atj  +  A<rij)njr(uf  +  Auf )dr, 

Jl'rr  Jr„ 

r  /-U^  +  L 

Jrcr  Ju\r-u f 


tfhd{uf-  -  uf)dTc 


(4.2) 
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where  .B  =  4<r  :  S  :  cr  is  the  complimentary  energy  density  and  AB(<7y,  Aery)  is  its  increment  due  to  stress 

l  2 

increase.  S  is  the  material  compliance  matrix.  The  notations  (•)  and  (•)  represent  two  sides  of  the  internal 
cohesive  crack  surface.  The  last  term  provides  the  work  done  by  the  cohesive  tractions  t(oh  due  to  crack 
surface  separation.  In  VCFE  formulation,  the  equilibrium  conditions  and  constitutive  relations  in  the  matrix 
and  the  compatibility  conditions  on  the  element  boundary  and  crack  surface  are  satisfied  a-priori  in  a  strong 
sense.  The  element  kinematic  equation: 

Vue  =  ee  in  (4.3) 

is  however  satisfied  in  a  weak  sense  from  the  stationary  condition  of  the  element  energy  functional  in 
equation(4.2).  The  weak  form  is  obtained  by  setting  the  first  variation  of  IIe  with  respect  to  stress  in¬ 
crements  to  zero,  i.e. 


-  f  (  +eij)  SAa,j  dfl+  [  SAcrij  ne-  ( uf  +  A uf)  dd(lc 

Jet e  o^ij  JdSie 

-  /,  SAoij  nf  («f  +  Auf)  dVcr  -  L  SAaij  n f  (uf  +  Auf )  dTcr  =  0  (4.4) 

Jrc  r  Jrcr 


Solution  of  equation  (4.4)  yields  domain  stresses.  Furthermore,  the  VCFE  formulation  assumes  weak  sat¬ 
isfaction  of  the  traction  reciprocity  conditions  on  (i)  the  inter-element  boundary  rme,  and  (iii)  the  domain 

1  2 

traction  boundary  rte  and  (iii)  the  crack  surfaces  Fcr  and  Ter: 


(ffij  +  A/Jij )n^+  —  —  ( <?ij  +  A(Xij)nf~  on  Tme  (inter-element  boundary) 

(aij  +  Aaij)nf  =  ti  +  Ati  on  rte  (traction  boundary) 

(aij  +  Aoij)lnf  =  (aij  -|-  Aaij)2nf  on  rcr 


(4.5) 


In  the  variational  principle,  the  weak  form  is  obtained  by  setting  the  first  variation  of  the  total  energy 

N  1  2 

functional  II  =  Y2e=i  with  respect  to  the  displacements  AuE,  Aucr  and  Aucr  respectively,  to  zero,  or 
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Oij  +  Aaij)n.j6uf  ddVl  - 


V  <5uf  e  vf  =  {vf  €  W°(Olf ) 


(<i  +  At*)  ]  <5uf  drtm  =  0 
=  0  on  rue}  V  e  on  dCle 


(4.6) 


and 


+  A  0ij)nf 

+  Aaij)rijr 


-  #]  Suf  dVcr  =  0 


+  4>'i\5<  dTcr  =  0 


V  (5Ugr  €  V£r,  V  e  on  Tcr 


(4.7) 


lr  ir-  lr  Ir  12 

where  rf>  =  +^u<  _u‘  _Au‘  t^ohd{u1r  -  uf)  is  the  cohesive  energy  function  and  c/>( 


..cr  ,,cr 


4.2.1  Cohesive  zone  models  for  crack  propagation 


dj, 

du\r  ' 


Cohesive  zone  models,  introduced  in  [6,  29]  and  developed  in  [71,  72,  73,  104,  34,  35,  42,  76],  are  effective  in 
depicting  material  failure  as  a  separation  process  across  an  extended  crack  tip  or  fracture  process  zone.  In 
these  models,  the  tractions  across  the  crack  reach  a  maximum,  subsequently  decrease  and  eventually  vanish 
with  increasing  separation  across  the  crack.  The  cohesive  model  used  in  this  chapter  is  a  three  parameter 
rate  independent  linear  cohesive  model,  proposed  in  [42,  76].  This  is  an  extrinsic  (two  stage)  model  which 
has  an  infinite  stiffness  or  slope  in  the  rising  portion  of  the  traction-separation  law  up  to  a  peak  traction 
value.  This  is  followed  by  linear  descending  segment  till  a  zero  traction  value  is  reached.  The  model  assumes 
a  free  cohesive  energy  potential  cj>  such  that  the  traction  across  the  cohesive  surface  is  expressed  as: 


t 


coh 


__  dcp 

ddnn  +  d6t 


(4.8) 
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Here  5n  and  5t  correspond  to  the  normal  and  tangential  components  of  the  opening  displacements  over  the 
cohesive  surface  in  the  n  and  t  directions  respectively.  An  effective  opening  displacement  is  defined  as 


6=y/%+fP6  ?  (4.9) 

where  (3  is  a  coupling  coefficient  to  allow  assignment  of  different  weights  to  normal  and  tangential  opening 
displacements.  Consequently  the  cohesive  surface  traction  reduces  to 

tcoh  —  ~(P26t  t  +  <5„n),  where  t  =  ~  =  J  f“h2  +  (3~2tfoh2  (4-10) 

o  oo  v 

where  t™h  and  tctnh  are  the  normal  and  tangential  components  of  surface  tractions.  The  effective  cohesive 
force  t  in  this  model  for  increasing  <5  takes  the  form 

V<5  <  8e  (4.11) 

0V<5  >  6e  (4.12) 

Se  corresponds  to  the  separation  at  which  t  goes  to  zero  and  <7mal  is  the  peak  value  of  t.  The  effective  normal 
traction-separation  response  of  this  model  is  depicted  in  figure  (4.2).  In  the  softening  region,  Unloading  from 
any  point  on  the  traction-separation  curve,  proceeds  along  a  linear  path  from  the  current  position  to  the 
origin  as  shown  by  the  line  BO  in  figure  4.2.  The  corresponding  t  —  5  relation  is 

-  Smax5  V<5  <  5max  <  5e  (4.13) 

Umax 

Reloading  follows  the  path  OBC  with  a  reduced  stiffness  in  comparison  with  the  original  stiffness.  Traction 
vanishes  for  S  >  Se. 

For  negative  normal  displacement  (compression),  stiff  penalty  springs  with  high  stiffness  are  introduced  be¬ 
tween  the  node-pairs  on  the  crack  face.  To  define  the  tangent  stiffness  matrix,  it  is  necessary  to  distinguish 
between  crack  initiation  (S  =  0)  and  crack  propagation  from  an  initialized  state  (<5  >  0).  In  the  former, 
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tc°h  =  t,  and  t1°h  =  0  are  assumed,  which  implies  that  the  initiation  is  in  pure  mode  I.  The  cohesive  pa¬ 
rameters  in  this  study  are  calibrated  from  experiments  done  for  epoxy-steel  composites  as  discussed  in  [58,  37]. 


Recent  experimental-computational  studies  on  composites,  conducted  in  [100]  show  that  the  three  or  four 
parameter  cohesive  models  are  more  suitable  for  modeling  interfacial  debonding  in  comparison  with  the  two 
parameter  models  based  on  Ferrante’s  law  [71,  72,  73].  Similar  conclusions  have  also  been  drawn  in  the  work 
by  Ghosh  et.  al.  [37,  58],  where  bilinear  cohesive  models  were  chosen  to  study  interfacial  debonding  in  fiber 
reinforced  composites. 


4.2.2  General  element  assumptions  and  weak  form 

In  the  absence  of  body  forces,  two  dimensional  stress  fields  satisfying  equilibrium  relations  can  be  generated 
from  the  Airy’s  stress  function  <h(a;,  y).  In  the  incremental  formulation,  stress  increments  are  obtained  from 
derivatives  of  the  stress  functions  A$(a :,y)  as: 


(  a  \ 


A(7v 


\A°*V  J 


(  d2A$  ^ 

a2A<i» 
dx 2 

d2A<E 

^  OxOy  J 


[P(x,y)]{A/3} 


(4.14) 


where  {A/3}  is  the  column  of  unknown  stress  increment  coefficients,  associated  with  the  stress  interpolation 
matrix  [P(x,j/)j.  Convergence  properties  and  efficiency  of  VCFEM  depend  on  the  choice  of  <F.  These 
functions  should  adequately  account  for  the  geometry  and  location  of  the  heterogeneity  in  the  element. 
Polynomial  functions  alone  do  not  contribute  to  this  requirement  and  hence  lead  to  poor  convergence  [68,  88]. 
Consequently,  stress  functions  in  X-VCFEM  are  constructed  from  different  expansion  functions  that  have 
complementary  effects  on  the  solution  convergence  for  the  propagating  crack.  Compatible  displacement  fields 
satisfying  inter-element  continuity  on  the  element  boundary  <9f if  and  intra-element  continuity  on  the  crack 
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face  Tcr  are  generated  by  interpolation  of  nodal  displacements,  [37,  68,  88]  as: 


{Aue}  =  [L,,]{Age}  on  dCle 

{A^CT}  =  [Lcr]{A<fr}  on  rcr, 

{Aucr}  =  [Lcr]{Aqcr}  on  T^  (4.15) 

1  2 

The  interpolation  matrices  [Le],  [L^],  [Lcr]  for  the  nodal  displacements  on  the  respective  boundaries  are 
constructed  using  standard  linear  or  hierarchical  shape  functions. 

Remark:  It  is  desirable  that  the  displacement  interpolations  on  the  crack  surface  in  equation  (4.15)  have  ad¬ 
equate  resolution,  consistent  with  the  high  resolution  in  the  stress  fields  near  the  crack  tip.  To  accommodate 
this,  hierarchical  shape  functions  are  added  to  standard  linear  shape  functions  to  describe  displacements  on 
the  crack  surface  as: 

4 

ucr  =  53  Ni(s)  *  qf  (4.16) 

»=i 

where  Ni  =  |(1  —  s),  N2  —  ^(1  +  s),  N$  =  i(s2  -  1),  and  N4  =  |(s3  -  s).  The  first  two  are  the  standard 
linear  shape  functions,  while  the  last  two  are  the  hierarchical  shape  functions  in  natural  coordinates  s.  The 
degrees  of  freedom  corresponding  to  higher  order  shape  functions  (i.e.  to  quadratic,  cubic,  etc.)  cannot  be 
interpreted  as  nodal  values  of  displacement.  Instead,  they  are  values  of  some  higher  order  derivatives  of  the 
solution  at  the  midpoints  (or  linear  combination  of  these  derivatives). 

Substituting  the  interpolations  of  stress  and  displacement  fields  from  equations  (4.14)  and  (4.15)  into 
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equation  (4.2)  results  in  the  matrix  form  of  the  element  complimentary  energy 


where 


He  =  -1-{f3  +  A(3}TlU]{(3  +  Ap}  +  {p  +  Af3}T{G}e{qe+Aqe} 


+  {P  +  A/3}T(Gcr]{qCr  +  Aqcr}  -  {t}T{qe  +  Aqe} 

-  {0  +  A/3}T[GCT']{q2cr  +  Aq^} 

112  2 
/*ucr+Aucr-ucr-Au 


/  /“  • 

Jrcr  Jucr-uc 


-  ir)dTc 


(4.17) 


[H]  = 


[Gc 


{t} 


/  [P]T[S][P]dfi  [Ge]  =  /  [P]T[ne][Le]ddfl 
J  n„  Janc 

(  !P]T(nCT][L1cr)rfrcr  [G-]  =  /  (P]T[ncr][Lcr]n!rc 
Jrcr  Jrcr 

[  {t  + At}r[Le]dr(m 
J  r(m 


(4.18) 


Construction  of  appropriate  stress  functions  with  optimally  high  resolution  is  necessary  for  accurately  de¬ 
picting  high  stress  gradients  near  the  crack  tip. 


4.2.3  Stability  conditions 

Following  the  stability  conditions  derived  for  displacement-based  and  stress-based  finite  element  approxima¬ 
tions  in  [5,  15,  111],  the  stability  conditions  of  the  stress-displacement  field  variational  problem  in  X-VCFEM 
depend  on  the  following  conditions. 

•  The  matrix  [H]  should  be  positive  definite.  From  the  definition  of  [H]  in  equation  (4.18),  the  necessary 
condition  for  it  to  be  positive  definite  is  that  the  compliance  tensor  [S]  be  positive  definite,  which  is 
true  for  elastic  problems. 

•  A  second  condition  is  that  the  finite-dimensional  stress  subspaces  T  be  spanned  uniquely  by  the  basis 
functions  [P] .  This  is  satisfied  by  assuring  linear  independence  of  the  columns  of  basis  functions  [P], 
which  also  guarantees  the  invertibility  of  [H], 
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•  Additional  stability  conditions  should  be  satisfied  to  guarantee  non-zero  stress  parameters  /3  for  all 
non-rigid  body  displacement  fields  on  the  element  boundary  uf  or  on  the  crack  face  uf.  This  is 
accomplished  by  careful  choice  of  the  dimensions  of  the  stress  and  displacement  subspaces,  i.e.  np  > 
nq  +n"  * 2  —  3,  where  np  is  the  number  of  /3  parameters,  and  np  and  n"  are  the  number  of  displacement 
degrees  of  freedom  on  the  element  boundary  and  crack  face  respectively. 

4.3  Creation  of  Enriched  Stress  Functions  in  X-VCFEM 

VCFEM  formulations  for  micromechanical  analysis  of  heterogeneous  materials  have  incorporated  polynomial 
and  reciprocal  stress  functions  based  on  analytical  micromechanics  results  in  [37,  68,  88,  38].  In  the  present 
work,  the  heterogeneity  is  in  the  form  of  an  evolving  cohesive  crack.  Two  conditions  need  to  be  considered 
in  the  choice  of  stress  functions.  The  first  is  that  it  should  adequately  represent  crack  tip  high  stress 
concentration  as  required  by  the  cohesive  zone  models.  Polynomial  functions  alone  are  unable  to  satisfy 
this  requirement  and  hence  suffers  from  poor  convergence.  The  second  condition  is  that  the  stress  function 
should  account  for  stress  jump  across  the  crack  surface.  The  stress  functions  in  X-VCFEM  incorporate  three 
different  components,  namely:  (a)  a  purely  polynomial  function  <f>po,y  to  yield  the  far  field  stress  distributions 
away  from  the  crack  tip,  (b)  a  branch  function  <t>b,'anch  that  is  constructed  from  level  set  functions,  and  (c) 

a  multi-resolution  wavelet  function  $wvlt  to  account  for  the  moving  crack  tip  stress  concentration.  Thus, 
$  =  <2>Poiy  4.  ^branch  <j ywvlt  _ 


4.3.1  Pure  Polynomial  Forms  of  Stress  Function: 

The  pure  polynomial  component  of  the  stress  function  $P0,y  is  written  in  terms  of  scaled  local  coordinates 
(£  =  * ~lc  i  V  ~  or'g>n  a!  the  element  centroid  { xc,yc ),  as: 

*po,y(^)=  P£  evqPpo  (4.19) 

p=0,q—0 

The  scaling  parameter  in  the  coordinate  representation  is  L  =  i/max(x  -  xc)  x  max(j/  -  yc) 

V(x,  y)  €  dQe.  The  use  of  the  scaled  local  coordinates  (£,  77),  as  opposed  to  global  coordinates  (x,  y)  in  the 
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construction  of  stress  functions,  prevents  ill  conditioning  of  the  [H]  matrix  due  to  the  high  exponents  of 
(x,y)  in  $P0,y.  As  discussed  in  [96],  invariance  of  stresses  with  respect  to  coordinate  transformations  can 
be  ensured  by  a  complete  polynomial  representation  of  $poiy,  while  stability  of  the  algorithm  requires  linear 
independence  of  the  columns  of  stresses  derived  from  <Fpoiy. 

4.3.2  Branch  Stress  Functions  Using  Level  Set  Methods 

The  branch  function  $>branch  facilitates  jumps  in  stresses  across  the  crack  surfaces.  These  functions  should 
not  affect  the  solutions  in  the  continuous  region  beyond  the  crack.  This  construction  requires  a  functional 
representation  of  the  surface  or  line  of  discontinuity.  Level  set  methods,  introduced  by  Sethian  [2,  93]  for 
following  the  evolution  of  interfaces,  is  ideal  for  representing  arbitrary  contours.  The  method  has  been  used 
by  Belytschko  and  coworkers  in  [12]  for  the  construction  of  branch  functions  associated  with  the  partition  of 
unity  in  a  displacement  based  FEM  formulation.  The  standard  level  set  methods  invoke  continuous  evolution 
of  the  entire  surface  of  discontinuity.  However  for  problems  involving  cracks,  the  only  evolution  occurs  at  the 
crack  tip  and  the  crack  surface  needs  to  be  frozen  behind  tip.  A  vector  level  set  method  has  been  developed  in 
[108,  107]  to  freeze  the  crack  surface  in  accordance  with  geometric  updating.  This  method  is  used  in  this  work. 

An  approximation  to  the  crack  surface  Tcr  in  figure  4.1  is  constructed  to  describe  the  discontinuous 
stress  fields  across  crack  paths.  As  shown  in  figure  4.3(a),  the  discontinuous  surface  is  expressed  by  a  signed 
distance  function  /(x)  defined  as 


/(x)  =  min  ||  x  -  x  ||  sign(n+  •  (x  —  x))  (4.20) 

where  x  is  a  point  on  the  surface  of  discontinuity  and  n+  is  a  unit  normal  pointing  in  the  direction  of  the 
region  of  positive  distance  function. 

Consequently,  x  is  the  closest  point  projection  of  any  point  x  on  Fcr.  In  order  to  describe  the  crack 
path  accurately,  the  signed  function  /(x)  is  evaluated  at  every  integration  point  in  the  Voronoi  cell  element 
directly.  The  process  of  constructing  branch  functions  involves  steps  that  are  described  below. 


61 


•  Radial  distance  functions  to  the  two  crack  tips  ri(x)  and  ^(x)  and  the  corresponding  angular  positions 
0i(x)  and  #2(x)  are  depicted  in  figure  4.3(a).  These  functions  are  expressed  in  terms  of  coordinates  of 
local  systems  (£,77)  with  origins  at  the  crack  tips.  For  the  local  system  at  crack  tip  1,  the  coordinates 
of  x  are  (£1,771).  In  accordance  with  the  definition  of  the  signed  distance  function,  the  radial  distance 
and  angle  functions  are  expressed  as 


n(x)  =  yJtf+Vi  and  0i  (x)  =  j 


,-1 L 


*  -  sin  £1  <  0,/  >  0 


—  Kin  * 


-sm  ~  7r  £i<0,/<0 


sin  £1  >  0 

rj  - 


(4.21) 


Similarly,  the  radial  distance  and  angle  functions  for  the  coordinate  system  at  crack  tip  2  are  defined 


as: 


r2(x)  =  \j (2  +  rti  and  02(x)  =  | 


■n  -  sin  1  £  £2  <  0,  /  >  0 

-sin~1^-n  &  <  0,  /  <  0 

fr>0 


(4.22) 


The  branched  stress  function  is  constructed  in  terms  of  the  functions  /(x),  0i,  n,  02,  and  r2,  as: 


^ branch 


s—0,t—0 


2  •  ^1  . 
r1stn^r. 


^COSY^lVlPst 


(4.23) 


The  terms  r\  and  r|  in  <pbranch  are  necessary  for  avoiding  crack  tip  singularity  in  the  stresses  due  to 
this  function  and  for  improving  the  accuracy.  Along  the  tangential  extension  to  the  crack  path  at  the 
tip  1,  $branch  is  zero  since  sin^-  =  0.  Hence  $i"'anch  does  not  contribute  to  the  stresses  ahead  of  the 
crack  tip  1.  In  an  analogous  manner,  $ branch  gQes  j-0  zero  a]ong  the  extension  to  the  crack  path  at 
the  tip  2,  since  cos^  =  0.  Therefore  $branch  does  not  contribute  to  the  stresses  in  this  region  also. 
However,  along  the  crack  surface  between  the  two  crack  tips,  sin^-  =  ±1  on  both  sides  of  the  crack,  and 
cos =  1.  This  renders  <$branch  in  equation  (4.23)  discontinuous  across  the  crack  path.  In  $branchf  Ql 
is  used  to  create  the  discontinuity  across  the  crack  surface,  while  02  eliminates  the  discontinuity  ahead 
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of  crack  tip  2.  In  some  special  instances  with  only  one  crack  tip  such  as  a  panel  with  an  edge  crack, 
equation  (4.23)  may  be  simplified  by  removing  r 2  and  02  dependence  to  yield 

^branch  =  £  r^sin^lVlPst  (4.24) 

S,< 


A  coordinate  transformation  is  required  to  obtain  stress  components  in  the  global  coordinate  system  from 
$ branch (^  ^)  based  on  the  local  coordinate  system. 


where  [Qi,]  is  the  transformation  matrix  from  (£1,771)  and  (£2,772)  systems  to  ( x,y ),  and  is  expressed  as 


din  di?i 
dy  dx 

d£i  as*. 

dy  dx 

dm  dm 

dy  dx 

d£ 2  d£ 2 
dy  dx 


(4.26) 


The  branch  function  is  evaluated  at  every  integration  point  in  the  element.  A  typical  function  <&branch 
for  s  =  0  and  t  =  0  is  plotted  in  figure  4.3(b).  The  plot  shows  that  the  function  is  continuous  everywhere  in 
the  domain  except  across  the  crack  surface.  The  example  of  a  double  cantilever  beam  under  a  sliding  load, 
as  shown  in  figure  4.7,  explains  the  effect  of  level-set  method  based  branch  functions.  In  figure  4.7(a),  the 
dimension  is  a  =  1.5m.  Figure  4.7(b)  shows  the  stress  axx  plots  as  a  function  of  y  at  x  —  —0.3 m.  The  stress 
functions  are  constructed  with  and  without  branch  functions  in  this  example.  axx  changes  its  sign  with  a 
jump  in  its  magnitude  on  different  sides  of  the  crack  and  the  jump  at  y  =  0  is  predicted  well.  However,  the 
transition  is  gradual  from  negative  to  positive  values  for  the  curve  without  branch  functions.  Although  the 
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transition  takes  place  in  a  short  interval,  the  method  is  not  able  to  catch  the  discontinuity  without  branch 

1  2 

functions.  This  also  results  in  the  matrices  [G"}  and  [GCT]  in  equation  (4.17),  on  different  sides  of  the  crack 
to  be  linearly  dependent  on  each  other  (one  is  the  negative  of  the  other). 

4.3.3  Multi-resolution  Wavelet  Functions  for  Modeling  Cohesive  Cracks 

Wavelet  bases,  discussed  in  [20,  70],  are  L2(1Z)  and  generally  have  compact  support.  Only  the  local  coeffi¬ 
cients  in  wavelet  approximations  are  affected  by  abrupt  changes  in  the  solution,  such  as  for  shock  waves.  This 
localization  property  makes  the  wavelet  basis  a  desirable  tool  for  problems  with  a  high  solution  gradients, 
concentrations  or  even  singularity.  A  brief  introduction  to  wavelet  basis  functions  is  provided  next. 

Principles  of  wavelets  and  multi-resolution  analysis 

The  construction  of  wavelet  functions  starts  from  a  scaling  or  dilatation  function  <p{x)  and  a  set  of  related 
coefficients  {p(k)}kez  which  satisfy  the  two-scale  relation 

#c)  =  £>(*00(2  x-k)  (4.27) 

k 

The  scaling  function  has  a  compact  support  only  if  many  coefficients  p(k)  are  non-zero.  Translations  of  the 
scaling  function  cf>(x  —  k )  form  an  unconditional  basis  of  a  subspace  Vo  C  L2(1Z).  Through  a  translation  of 

4>  by  a  factor  of  2"  and  dilation  by  a  factor  of  k  ■  2~n  the  unconditional  basis  is  obtained  for  the  subspace 

Vn  C  L2{Tl)  as 

<t>n,k(x)  =  2n'2<j>(2nx-k)  (4.28) 

for  a  resolution  level  n.  The  scaling  function  <j>  is  defined  as  orthonormal  if  translations  at  the  same  level  of 
resolution  satisfies  the  condition 


/OO 

<t>n,k(x)<l)n,i(x)dx  =  5k,i  V  n,k,l  e  z 

-OO 


(4.29) 
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Consequently,  the  best  approximation  of  a  function  f(x )  in  the  subspace  Vn  of  L2( 72.)  is  expressed  as  the 
orthogonal  projection  of  /  on  Vn  as: 


/OO 

f(x)<t>n,k(x)dx  (4.30) 

K  ■ °° 

Approximation  of  f(x),  can  be  made  at  different  resolution  levels,  and  these  approximations  in  subspaces 
•••,  V„-i,  Kj,  Cu+i,  follow  the  relation 

{0}  =  VLoo  C  •  •  •  C  V_!  C  V&  C  Vi  C  •  •  •  C  Voo  =  L\n),  where 

lirrin^ooVn  =  (^J  Vn  is  dense  in  L2(TZ)  and  fimn__0O  fl„  V„  =  {0}  (4.31) 

In  the  multi-resolution  level  transition,  the  information  lost  in  the  transition  from  level  Vn+\  to  level  Vn  is 
characterized  by  an  orthogonal  complementary  subspace  Wn.  A  basis  for  the  subspace  Wn  can  be  obtained 
is  in  the  same  manner  as  for  scaling  function,  i.e.  by  dilating  and  translating  the  mother  wavelet  function 

V-’(z)  =  <?(fc)^(2a;  -  k)  (4.32) 

k 

The  subspaces  spanned  by  the  wavelet  functions  have  the  following  essential  properties: 

(i)  Vn+ 1  =  Vn  ©  Wn  V,  i.e.  Wn  is  the  orthogonal  complement  of  Vn  toV^+i; 

(«)  For  orthonormal  bases,  Wn\  is  orthogonal  to  Wn2\ 

(in)  For  orthonormal  bases,  3)'^L_00  Wn  =  L2(1Z)  (4.33) 


An  approximation  of  the  function  f(x)  at  the  n  —  th  resolution  level  may  be  expressed  as  the  orthogonal 
projection  of  /  on  Wn  as 


/  ->  Vnf  (x)  = 

k 


/oo 

f(x)lpn:k(x)dx 

-oo 


(4.34) 
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Due  to  the  orthonormality  and  multi-resolution  properties  of  wavelet  basis  functions,  higher  level  approximate 
solutions  can  be  generated  from  results  of  lower  level  solutions  (see  [20,  70])  by  selective  superposition  of 
complementary  solutions.  The  use  of  adaptive  enrichment  is  very  attractive  to  those  regions  where  a  pre¬ 
determined  ’error  or  residual’  tolerance  is  not  met  at  the  lower  level. 

Selection  of  the  wavelet  function 

Various  wavelet  functions  have  been  proposed  in  the  literature  for  numerical  solutions  of  ODEs  and  PDEs. 
These  functions  have  been  incorporated  in  the  method  of  weighted  residuals  like  the  Galerkin’s  method  and 
collocation  method  to  solve  problems  with  multi-level  features  in  [41,  51,  80],  Among  the  large  number  of 
wavelet  functions  proposed  are  the  Haar  function  [43],  the  Meyer’s  wavelets  [65],  the  Chui- Wang’s  B-spline 
wavelets  [21],  etc.  One  of  the  most  commonly  used  wavelet  functions  is  Daubechies’  compactly  supported 
orthonormal  wavelets  [25,  27,  41].  However,  they  are  constructed  through  recursive  algorithms  and  do  not 
have  an  explicit  analytic  expressions.  This  makes  it  is  difficult  to  obtain  their  first  and  second  derivatives, 
which  is  a  requirement  in  X-VCFEM  for  deriving  stresses  in  terms  of  stress  functions.  Also  the  orthonormality 
of  the  Daubechies  wavelet  cannot  be  transferred  to  the  orthonormality  for  stresses  by  differentiation,  and 
hence  they  are  not  considered  to  be  suitable  for  stress  functions  in  X-VCFEM.  Alternatively  a  family  of 
Gaussian  functions,  for  which  the  first  and  second  order  derivatives  are  popular  wavelets  bases  [14,  31,  57], 
is  implemented  in  the  representation  of  X-VCFEM  stress  functions  and  stresses.  The  expressions  for  the 
Gaussian  function  and  its  n-th  order  derivative  are: 


G{x)  —  e  fa)/2  and  =  (—  l)n-^~ (e  Co6)2/2)  (4.35) 

The  dilation  and  translation  parameters  a  and  6  respectively  can  assume  arbitrary  values  and  can  be  changed 
in  a  continuous  fashion.  The  ability  of  wavelets  to  translate  diminishes  the  need  to  re-define  new  elements 
or  remesh  in  conventional  FEM  solution  of  problems  with  moving  boundaries.  By  changing  translation 
parameters,  the  multi-levels  of  wavelet  bases  can  be  made  to  closely  follow  a  moving  crack  tip.  Additionally 
the  dilation  parameter  with  compact  adjustable  window  support  can  be  used  to  provide  high  refinement 
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and  resolution.  Hence  it  is  a  convenient  way  of  moving  the  stress  concentrations  using  the  multi-resolution 
properties. 

Multiresolution  wavelet  based  stress  functions  for  crack  problems 

The  wavelet  based  stress  function  is  constructed  in  a  local  orthogonal  coordinate  system  (£,??),  centered  at 
the  crack  tip.  The  £  direction  corresponds  to  the  local  tangent  to  the  crack  surface.  The  corresponding 
stress  function  $a,b,c,d  in  the  Gaussian  wavelet  basis  is  given  as: 


$a,b,cAZ<v)  =  e 


')V2&,6,C,d 


(4.36) 


where  a,b,c,d  are  parameters  that  can  take  arbitrary  continuous  values.  For  implementation  in  multi¬ 
resolution  analysis  involving  discrete  levels,  the  translation  and  dilation  parameters  should  be  expressed  as 
discrete  multiples  of  some  starting  values.  Consequently,  these  discrete  values  am,  bn,  Ck  and  d;  are  expressed 
as: 

dm  —  0,1  ■ 
bn  —  71 '  b\  ■  am 

<  (4.37) 

Ck  =  ci  ■  (trc)fc_1 

d\  —  l  •  d  i  •  Ck 

Here  ( m,k )  correspond  to  the  levels  and  (n,l)  correspond  to  the  discrete  translation  of  the  bases  in  the 
(G7?)  directions  respectively.  The  parameters  (ai,ci)  are  the  initial  dilating  values  at  the  first  level  m  —  1, 
while  tra(<  1),  trc(<  1)  are  the  transfer  rates  from  one  level  to  the  next  higher  one.  The  parameters  b\,d\ 
represent  the  starting  values  of  a  step  translation  quantity  at  the  m  —  th  dilation  level.  The  narrow  (higher 
level)  wavelets  are  translated  by  small  steps,  whereas  the  wider  (lower  level)  wavelets  are  translated  by  large 
steps.  Parameters  tra  =  trc  =  1  and  6i  =  d\  =  0  imply  no  dilation  and  translation  respectively.  Parameters 
cq,  cc,  and  do  are  counterparts  of  do,  ac,  and  bo  in  ?;  direction.  With  the  specific  relations  between  dilation 
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and  translation  parameters,  the  Gaussian  wavelet  enriched  stress  function  in  equation  (4.36)  becomes 


$ 


m,n,k 


,z  (£.*?) 


e-<^)V2e-( 


ck 


f/2 


m,ntktl 


(4.38) 


The  family  of  wavelet  enriched  stress  functions  in  equation  (4.38)  are  not  orthonormal,  but  they  construct 


a  linearly  independent  basis  [26].  This  leads  to  robustness  and  high  precision  in  the  reconstruction  of  any 


function  /  even  with  low  level  coefficients.  The  wavelet  enriched  stress  function  in  X-VCFEM  is  thus  written 


as 


2^  ,fcn  ,1-n 

Tn=l,n=-^-,fc=l,Z=0 


(4.39) 


The  corresponding  stresses  are: 
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(4.40) 


The  stress  components  in  the  global  coordinate  system  are  obtained  by  the  transformation  from  the  local 


coordinate  system  as 
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(4.41) 


where  [Qtu]  is  the  transformation  matrix  from  (£,  r])  to  (x,y): 


Figure  4.4  shows  the  support  region  for  the  wavelets  enriched  ^wvlt(^,r])  in  a  X-VC  element.  This  region  is 
positioned  symmetrically  in  the  vicinity  of  evolving  crack  tips.  The  crosses  (x)  corresponds  to  the  position 
of  each  wavelet  basis  function  bn,dn  at  a  lower  level,  and  the  squares  (□)  correspond  to  additional  locations 
at  a  higher  level  in  the  multi-resolution  algorithm.  Only  the  points  at  the  top  half  are  shown  in  the  figure 
due  to  symmetry. 

The  method  of  implementation  of  the  multi-resolution  wavelet  enriched  stress  functions  in  X-VCFEM  is 
described  below. 

1.  For  the  starting  level  m  =  k  =  1,  20  points  marked  by  crosses  (x)  in  figure  4.4  (a),  are  used  to  delineate 
the  wavelet  enriched  function  $wvlt(£,T])  in  equation  (4.39).  This  corresponds  to  m  =  1,  n  —  5,  k  =  1 
and  l  =  4. 

2.  With  ensuing  higher  levels  in  the  multi-resolution  wavelet  functions  according  to  the  equation  (4.37), 
higher  level  wavelet  bases  are  added  to  the  stress  function  as  marked  by  squares  (□)  in  figure  4.4 
(b).  The  addition  is  done  adaptively  in  accordance  with  error  criteria  discussed  in  section  (4.3.4).  A 
refinement  in  the  starting  region  of  wavelet  enrichment  occurs  in  each  added  level,  i.e.  the  window  size 
of  additional  wavelet  basis  functions  is  smaller  than  ones  at  a  lower  level.  This  allows  a  zoom  in  to 
catch  higher  gradients  that  are  missed  at  the  coarser  scales. 

3.  The  process  of  successive  multi-level  refinement  can  continue  till  a  predetermined  error  tolerance  is 
reached. 

Remark:  The  line  of  the  cohesive  crack  is  likely  to  intersect  the  region  of  support  of  the  wavelet  bases  func¬ 
tions.  It  is  important  for  the  numerical  algorithms  to  assure  that  wavelet  functions  based  on  one  side  of  the 
cohesive  crack  does  not  contribute  to  stresses  on  the  other  side.  The  influence  of  wavelet  stress  functions 
should  be  cut  off  across  this  line  of  discontinuity  by  establishing  a  truncated  effective  support  domain  for  the 
wavelet  function.  This  is  accommodated  by  ignoring  the  contribution  of  quadrature  points  in  the  numerical 
integration  on  the  other  side  of  the  crack  as  detailed  in  section  4.5.3. 
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In  summary,  the  stresses  in  an  element  are  computed  by  adding  contributions  from  equations  (4.21), 


(4.23)  and  (4.40),  to  yield 
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(4.43) 


4.3.4  Error  measure  for  adaptive  wavelet  enrichment 

The  Euler  equation  (4.3)  indicates  that  the  error  in  the  kinematic  equation,  which  is  satisfied  in  a  weak 
sense,  may  be  primarily  attributed  to  the  lack  of  adequate  resolution  in  the  equilibrated  stress  fields.  A 
strain  energy  based  element  error  measure,  derived  in  [88],  is  extended  to  the  present  problem.  Let  a  stress 
field  be  enriched  from  a  level  n  to  level  n  +  1  by  adding  the  wavelet-based  enrichment  stress  crenT ,  i.e. 


0.level(n+ 1)  __  ^level(n)  _|_  ^.enr 


(4.44) 


The  corresponding  percentage  change  in  the  strain  energy  ( SE  =  fn  OijSi3ki&ki  dQ),  may  be  expressed  as 


A  SE  = 


SE(<rUvel:n+1)  -  SE(crlevel:n) 


SE(cr 


level:n-\- 1 


x!00% 


(4.45) 


In  view  of  the  local  properties  of  wavelets  and  stress  concentration  at  crack  tips,  the  strain  energy  in 
equation  (4.45)  is  calculated  only  in  a  small  region  around  crack  tip  fienr-  Adding  levels  is  conditioned  upon 
the  requirement  that  A  SE  is  less  than  a  preset  tolerance,  which  in  this  work  is  chosen  to  be  «  4%. 
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4.4  Solution  Method 


Crack  growth  in  multiply  cracked  materials  is  solved  using  an  incremental  approach,  where  a  set  of  elemental 
and  global  equations  are  solved  in  each  increment  for  stresses  and  displacements. 


1.  Local  equations  for  each  element  are  obtained  by  substituting  the  stress  interpolations  of  equation  (4.43) 
and  boundary/crack  face  displacement  interpolations  of  equation  (4.15)  in  the  element  energy  functional 
equation  (4.17)  and  setting  its  variation  with  respect  to  the  stress  coefficients  A/3  to  zero.  This  results  in 
the  weak  form  of  the  element  kinematic  relations 


[H]e{/3  •+•  A/3}e 


[Ge]  [Gc 


-  [Gc 


qe  +  Aqe 

1  l 

qcr  +  Aqcr 

2  2 
qCT  +  Aqcr 


(4.46) 


or  in  a  condensed  form 


[H]e{/3  +  A/3}e  =  [Gje{q  +  Aq}e 


(4.47) 


Since  equation  (4.47)  is  linear,  the  stress  coefficients  can  be  directly  expressed  in  terms  of  the  nodal  displace¬ 
ments,  provided  the  element  [H]e  matrix  is  invertible. 


2.  Subsequently,  the  weak  forms  of  the  global  traction  continuity  conditions  are  solved  by  setting  the 
variation  of  the  total  domain  energy  functional  with  respect  to  the  generalized  displacement  components  to 
zero.  This  results  in  the  weak  form  of  the  traction  reciprocity  conditions 
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or  in  a  condensed  form: 


N 


N 
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(4.49) 
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The  forces  at  the  crack  surface  are  expressed  in  terms  of  the  cohesive  energy  as 
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(4.50) 


Combining  equations  (4.47)  and  (4.49)  and  eliminating  the  stress  coefficients  {/3  +  A /3}e,  results  in  the 
equation  for  solving  the  generalized  displacements 
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(4.51) 
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e=l 


Equation  (4.51)  is  a  nonlinear  matrix  equation  system  due  to  the  cohesive  laws.  Consequently,  a  Newton- 
Raphson  iterative  solver  is  invoked  to  solve  for  the  increments  of  nodal  displacements.  The  linearized  form 
of  equation  (4.51)  for  the  j-th  iteration  is 
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which,  in  a  condensed  form  is 


IKW  =  {RLJ  -  {Rfn 


(4.53) 


A  numerical  problem  associated  with  modeling  cohesive  crack  growth  is  the  occurrence  of  snap-back  as  is 
shown  in  the  macroscopic  load-deformation  behavior  plot  of  figure  4.5. 

This  has  been  discussed  for  a  three  point  bending  solution  in  [66].  For  a  deformation  controlled  process 
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with  monotonically  changing  deformation,  the  solution  ignores  the  reverse  portion  of  the  displacement  BCD, 
occurring  with  snap-back.  The  Newton-Raphson  solver,  where  the  loading  process  is  monotonically  controlled 
by  incremental  deformation  or  load  conditions,  exhibits  a  discontinuous  drop  from  point  B  to  point  D.  It 
is  obvious,  that  this  solver  needs  to  be  augmented  with  the  capability  to  account  for  the  part  BCD,  i.e.  to 
decrease  both  load  and  deformation  with  the  growth  and  opening  of  the  crack.  The  arc-length  solver  has 
been  proposed  in  [23,  24,  92]  as  a  method  of  overcoming  this  shortcoming  by  introducing  an  unknown  loading 
parameter  (A  +  d\)  to  govern  the  load  increments.  Equation  (4.53)  is  modified  with  this  loading  parameter 
as 


[KSW  =  (AJ  +  dXj){R9cxt}  -  {Rfnty  (4.54) 

where  both  4AJ  and  dq1  are  unknowns,  and  d\i  can  be  either  positive  or  negative.  The  additional  unknown 
dAJ  requires  the  solution  of  a  constraint  equation,  written  in  terms  of  the  magnitude  of  the  deformation  of 
all  the  nodes  on  the  crack  surface  as 


53  ((Auf)2  +  (Auf )2)  =  Ai2  (4.55) 

i£Crk 

where  Crk  represents  the  set  of  all  nodes  on  crack  surfaces.  A  summary  of  the  solution  process  is  explained 
in  the  flowchart  of  figure  4.6. 

4.5  Aspects  of  Numerical  Implementation 

4.5,1  Adaptive  criteria  for  cohesive  crack  growth 

A.  Direction  of  incremental  cohesive  crack  advance:  In  linear  elastic  fracture  mechanics,  it  is  common  to  use 
the  “maximum  hoop  stress  criterion”  to  determine  the  direction  of  crack  propagation  [9,  12].  Cracks  are 
assumed  to  propagate  in  a  direction  normal  to  the  maximum  hoop  stress  in  this  criterion.  Since  stresses  at 
crack  tip  are  singular  in  LEFM,  stress  intensity  factors  are  usually  used  to  determine  the  direction  of  crack 
propagation.  This  criterion  is  only  suitable  for  K-dominated  problems,  where  the  size  of  the  fracture  process 
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zone  is  small  compared  to  the  size  of  the  specimen.  A  different  criterion,  based  on  the  cohesive  energy  at 
the  crack  tip  is  used  in  X-VCFEM.  A  relation  between  the  cohesive  energy  <p  for  complete  decohesion  and 
the  critical  energy  release  rate  Gc  has  been  established  in  [76]  from  the  definition  of  the  J— integral  as: 


Gc 


j  fR  85  , 

=  J=  t  -K—axi  = 

Jo  oil 


tdS  = 


(4.56) 


where  R  is  the  length  of  the  cohesive  zone.  Consequently,  the  crack  growth  direction  is  estimated  as  that, 
along  which  Gc  or  equivalently  the  cohesive  energy  <j>  is  maximized  for  a  given  crack  tip  state  of  stress.  The 
cohesive  energy  4>a  at  the  crack  tip  A  along  any  direction  a  can  be  expressed  for  an  arbitrary  separation 
<5(a)  as: 

/  a(Q)  \  /  ftM  85  \ 

t(a)d8j  =  t(a)  ■  —dtj  (4.57) 

where  t{a)  =  \/ (t“/l)2  +  @~2(tc°h)2  is  the  magnitudes  of  the  effective  cohesive  force.  The  corresponding 
unit  normal  n  and  tangential  t  vectors  along  the  direction  a  are  expressed  as 


n  =  —sinai  +  cosaj  ,  t  =  cosai  +  sinaj 


(4.58) 


The  normal  and  tangential  components  of  the  cohesive  traction  force  at  an  angle  a  may  then  be  deduced  as: 
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(4.59) 


and  hence  the  effective  cohesive  traction  for  direction  a  is 


f(a)  = 


(<7xxsin2a  —  axysin2a  +  ayycos2a)2  +  0~2(—  -crxxsin2a  +  uxycos2a  +  -<7yysin2a)2 


(4.60) 


The  incremental  direction  of  crack  propagation  is  assumed  as  that  which  maximizes  the  cohesive  energy  at 
A,  according  to  the  criterion: 


dMa)  =  0  ,„d  30021  <  0 
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(4.61) 


A  combination  of  equations  (4.57),  (4.60)  and  (4.61),  yield 
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(4.62) 


Equation  (4.62)b  results  from  the  fact  that  t  cannot  be  equal  to  zero  for  decohesion  to  initiate  and  hence 
the  necessary  condition  evolves  from  its  derivative.  The  direction  of  crack  propagation  ac  is  obtained  as  the 
solution  of  equation  (4.62)b  as 
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The  optimal  angle  a*VCFBM  is  chosen  as  the  one  that  satisfies  the  condition  in  equation  (4.62)c.  The 
corresponding  angle  given  by  the  maximum  hoop  stress  criterion  in  LEFM  is  expressed  in  terms  of  the  stress 
intensity  factors  K[,  Ku  as: 
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LEFM 


=  2  arctan^  [k^Ku  ±  ^(Kt/Kn)2'^) 


(4.64) 


where  the  sign  is  chosen  to  make  the  hoop  stress  positive.  The  first  of  equation  (4.63),  which  is  the  only 
choice  for  0  —  1,  exactly  matches  the  angle  given  by  the  maximum  hoop  stress  criterion  (4.64).  For  the  pure 
sliding  problem  shown  in  figure  4.7,  ac  predicted  by  the  equation  (4.64)  is  70.5°,  while  that  by  X-VCFEM 
for  cohesive  stresses  is  68.2°. 


B.  Length  of  the  incremental  cohesive  crack  advance:  Upon  establishing  the  direction  of  incremental  cohe¬ 
sive  crack  growth  ac,  the  length  of  cohesive  zone  advance  (A I)  should  be  estimated  in  the  crack  evolution 
scheme.  The  criterion  used  is  that  the  cohesive  energy  goes  zero  at  the  end  of  the  new  segment  as  shown  in 
figure  4.8(a).  To  achieve  this,  the  cohesive  energy  at  two  points  A  (present  crack  tip)  and  B  (close  to  A  in 
the  direction  of  crack  propagation)  are  evaluated  by  substituting  the  stresses  in  equation  (4.62)a.  The  tip  of 
the  cohesive  zone  is  obtained  from  the  linear  extrapolation  of  this  line  to  yield  zero  cohesive  energy.  From 
figure  4.8  (a),  the  increment  of  cohesive  crack  length  is  defined  as: 
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4>a  -  <t>B 
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A  l  = 


(4.65) 


C.  Cracks  crossing  the  interelement  boundaries  and  merging  with  each  other:  Crack  advance  from  one  Voronoi 
cell  element  to  the  next  is  conducted  in  X-VCFEM  using  an  algorithm  depicted  in  figure  4.8(b).  A  contin¬ 
uous  tracking  method  is  implemented  to  monitor  if  a  cohesive  surface  has  reached  or  gone  past  an  element 
boundary.  In  this  method,  the  intersection  of  the  crack  surface  and  an  element  boundary  is  obtained  by 
solving  the  equation  system 

x  -  Xj  y-Vi  x-xn  _  y  -  yn 

Xj  1  X{  1  “  Vi  3-n+l  Xn  Vn  + 1  Vn 

where  ( Xi,yi )  represents  the  tip  of  the  cohesive  crack  line  for  the  ith  increment,  and  (xn,  yn)  is  the  position 
of  the  nth  node  on  the  element  boundary.  If  the  intersection  point  is  outside  of  the  cohesive  line  or  the 
element  boundary,  no  intersection  is  assumed.  Once  a  cohesive  crack  has  reached  its  intersection  with  the 
boundary,  a  new  node  pair  (n i  712)  is  introduced  on  the  element  boundary  at  this  point.  The  node  pair 
belongs  to  the  intersection  of  the  element  boundary  and  the  cohesive  crack,  i.e.  nin.2  €  dfif  PI  Tcr.  The 
crack  is  subsequently  advanced  to  the  next  element  following  the  usual  procedure  outlined  before. 

Another  condition  that  is  considered  in  this  work  is  the  merger  of  multiple  cracks  as  shown  in  figure  4.8(c) 
for  two  cohesive  cracks.  The  algorithm  for  crack  merging  is  an  extension  of  the  intersection  algorithm, 
discussed  above.  At  the  end  of  every  increment,  all  the  cracks  that  have  propagated  in  that  increment  are 
recorded.  Subsequently,  the  intersection  of  the  last  incremental  segment  of  the  cohesive  crack  with  those  of  all 
neighboring  cracks,  that  belong  to  either  the  same  element  or  neighboring  elements,  is  checked  using  equation 
(4.66).  Once  the  intersection  of  two  crack  segments  is  ascertained,  a  three-node  junction  (mi, m2, m3),  as 
shown  in  figure  4.8(c),  is  inserted  at  the  point  of  intersection.  The  contribution  of  the  junctions  nodes  e.g. 
(mi,  m2)  to  the  load  vector  in  the  assembled  matrix  equation,  requires  special  treatment.  For  each  of  these 
nodes,  contributions  of  integrals  from  adjoining  crack  segments  belonging  to  two  different  cohesive  cracks, 
are  summed. 


(4.66) 
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4.5.2  Evaluation  of  stress  intensity  factors 


The  stress  intensity  factors  and  J—  integral  are  evaluated  in  the  post-processing  phase  of  the  computations. 
From  linear  fracture  mechanics,  the  relation  between  J—  integral,  stresses  and  stress  intensity  factors  are 
given  as 


J  ■■ 


IA 


Gik^ik^lj  0’ijUi1\)7ljds  = 


Kh 

E * 


(4.67) 


where  E*  —  E  (Young’s  modulus)  for  plane  stress,  E*  —  jzA  f°r  plane  strain,  and  u  is  the  Poisson’s 
ratio.  In  displacement  based  FEM  [66,  30],  the  contour  integral  is  converted  into  a  domain  integral  to 
improve  the  accuracy  of  the  stress  intensity  factors,  since  the  stresses  are  more  accurate  in  the  interior  of  an 
element.  However  in  X-VCFEM,  stresses  on  the  contour  and  the  interior  are  equally  accurate  due  to  stress 
interpolation  and  the  contour  integral  can  provide  similar  accuracy  as  the  domain  integral.  A  method  to 
extract  the  stress  intensity  factors  Ki  and  Ku  from  the  J-integral,  proposed  in  Yau  [114],  is  implemented 
in  X-VCFEM.  Displacement  fields  are  not  interpolated  in  the  interior  of  the  Voronoi  cell  element,  and  hence 
the  term  1/2,1  in  equation  (4.67)  requires  a  special  evaluation  method. 

1.  Compute  €n,  £22,  and  £12  at  a  series  of  points  (xj,i/i,  i  =  in  a  small  shadowed  region  around 

the  integration  point  (xo,2/o)  in  figure  4.8(d).  The  displacement  gradient  ui,i  is  calculated  from  en. 

2.  For  evaluating  1/2,1  displacements  u\  and  U2  at  any  point  (xj,yt)  are  interpolated  using  polynomial 
functions, 


ui{xi,yi)  =  a0+aiXi  +  a2yi+a3Xi , 

u2{xi,yt)  =  b0  +  hxi  +  b2yi  +  b3x^  +  ■  ■  ■  (4.68) 

where  ao,  ai,  •  •  • ,  and  bo,  bi,  •  •  • ,  6m  are  unknown  coefficients.  To  constrain  the  rigid  body  motion, 
coefficients  are  evaluated  from  displacement  values  at  two  points  on  boundaries. 

3.  Displacement  gradient  expressions,  Ui,i  {xi,  yi),  112,2(^1,1/;),  ^i,2(xi,Pi)  and  t^.iOei,  2/i)  are  obtained  by 
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taking  derivatives  of  the  expressions  in  equation  (4.68).  Strain  expressions  in  terms  of  the  unknown 
coefficients  are  computed  from  these  derivatives.  At  each  point  (re*,  yi),  the  strains  can  also  be  computed 
from  the  known  stresses  and  the  compliance  tensor,  i.e.  {e}  =  [S] [P]  {/3}.  The  unknown  coefficients  ao, 
ai,  ■  ■  ■ ,  cim  and  b0 ,  blt  •••,6m  in  equation  (4.68)  are  estimated  by  solving  a  least  square  minimization 
problem  for  the  strains.  Subsequently  the  displacement  gradient  U2,i  is  determined  at  the  integration 
point  (xo,y0). 


4.5.3  Numerical  integration  schemes  for  [H]  and  [G]  matrices 

Integration  of  [H]  matrix: 

Numerical  integration  over  each  element  is  conducted  by  the  Gaussian  quadrature  method  to  form  the  matrix 
[H]  in  equation  (4.18).  In  this  method,  each  Voronoi  cell  element  is  recursively  subdivided  into  triangular 
subdomains,  on  which,  integration  points  are  generated  for  the  Gaussian  quadrature.  The  steps  involved  are 
discussed  below. 


1.  For  each  Voronoi  cell  element  shown  in  figure  4.9,  the  centroid  O  is  first  generated.  The  first  set  of 
triangular  subdomains  is  created  by  joining  each  of  the  vertices  of  the  cell  e.g.  (A,  B,  C,  D,  E, 
F)  with  the  centroid  O. 

2.  Each  triangle  is  further  subdivided  into  two  triangles  if: 


Area  of  triangle 
Area  of  Voronoi  cell  element 


>  TOLarea 


(4.69) 


For  the  subdomain  triangle  BCO  shown  in  figure  4.9,  two  triangles  are  created  by  bisecting  the  longest 
edge  BC  at  O’  and  joining  it  with  the  opposite  vertex  O. These  new  smaller  triangles  are  again  checked 
against  the  tolerance  condition  and  further  dissection  is  executed  if  necessary.  Numerical  integration 
in  each  triangular  subdomain  is  done  using  13  Gauss  points. 

3.  For  the  region  containing  the  crack  tip  shown  in  figure  4.9,  a  smaller  value  of  TOLarea  is  chosen  in 
comparison  with  other  regions.  This  facilitates  a  higher  density  of  integration  points  in  regions  of  high 
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stress  gradients.  The  tolerance  in  an  element  is  consequently  adjusted  according  to  the  distance  of  the 
center  of  the  triangular  subdomain  from  the  crack  tip,  i.e. 


TOL 


area 


=  TOL 


min 

area 
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(TOL 


max 

area 


-  TOL™)  *  dtri 
L 


(4.70) 


where  L  is  a  scaling  parameter  defined  in  subsection  (4.3.1),  dtTl  is  the  distance  of  the  crack  tip  from 
the  subdomain  and  TOL™rae*,TOL™™a  are  assumed  tolerances.  In  this  work  the  tolerances  are  chosen 
as  TOL™xa  =  10%  and  TOL™"  =  1%. 

4.  The  intersection  of  the  support  of  wavelet  functions  with  the  cohesive  crack  line  call  for  a  truncated 
support.  This  is  done  by  eliminating  the  contribution  of  quadrature  points  that  lie  on  the  other  side 
of  crack  face  from  the  wavelet  center.  A  visibility  criterion  introduced  in  [10]  provides  an  easy  way  to 
accommodate  this  discontinuity  in  the  construction  of  truncated  support.  In  this  method,  the  cracks 
are  considered  to  be  opaque  when  generating  valid  numerical  integration  regions.  A  ray  is  emitted 
from  the  center  W  of  a  wavelet  basis  function  in  an  arbitrary  direction  as  shown  in  figure  4.4.  If  it 
encounters  an  internal  crack,  the  ray  is  terminated.  All  quadrature  points  lying  in  the  dark  shadow 
region  on  the  other  side  of  the  crack  CC'  are  suppressed  during  numerical  integration  of  this  wavelet 
basis. 


Integration  of  the  [G]  matrices: 

1  2 

In  equation  (4.18),  the  matrices  [Gcr]  and  [Gcr]  are  numerically  integrated  over  the  crack  surfaces  and  the 
matrix  [Ge]  over  the  element  boundary.  All  numerical  integrations  on  the  element  boundary  and  crack 
surfaces  are  executed  using  the  Gaussian  quadrature  method.  The  number  of  integration  points  Nint  on 
each  boundary/crack-face  segment  depends  on  the  distance  ds ide  between  its  center  and  the  crack  tip,  and 
is  chosen  from  the  condition 


Nint  =  < 


9 


dside  ^  0.1L 


16  dside  <0.1  L 


(4.71) 
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where  L  is  the  scaling  parameter. 


4.5.4  Invertibility  of  the  [H]  matrix 

A  nonsingular  or  invertible  [H]  matrix  necessitates  the  linear  independence  of  the  columns  of  the  [P]  matrix. 
For  pure  polynomial  expansions  of  the  stress  functions,  this  condition  is  naturally  attained.  However  when 
adding  the  other  terms,  some  of  the  terms  in  the  branch  and  wavelet  functions  may  have  linear  dependence 
on  the  polynomial  terms.  In  X-VCFEM,  the  rank  of  the  [P]  matrix  is  first  determined  from  the  diagonal 
matrix  resulting  from  a  Cholesky  factorization  of  the  square  matrix 

[H*j  =  /  [P]T[P]dn  (4.72) 

Nearly  dependent  columns  of  [P]  will  result  in  very  small  pivots  during  Cholesky  factorization.  The  cor¬ 
responding  branch  and  wavelet  function  terms  are  dropped  from  the  stress  function  to  prevent  numerical 
inaccuracies  in  inverting  (H) . 


4.5.5  Elimination  of  element  rigid  body  modes 

X-VCFEM  uses  a  stress-based  formulation  with  independent  representation  of  displacement  fields  on  the 

element  and  crack  boundaries.  In  general,  the  nodes  of  the  crack  face  are  not  topologically  connected  to 

the  element  boundary  nodes.  However  it  is  important  that  all  nodes  in  the  element  possess  the  same  rigid 

body  modes.  The  rigid  body  modes  of  the  element  boundary  displacements  {qe}  are  directly  constrained 

in  the  solution  process  through  prescribed  displacement  boundary  conditions.  However,  it  is  necessary  to 

1  2 

connect  these  with  rigid-body  modes  for  the  crack  face  displacement  fields  {qcr}  and  {qcr}.  Singular  value 
decomposition  or  SVD  has  been  discussed  in  [88]  as  an  effective  method  for  identifying  and  constraining 
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rigid  body  modes  at  interfaces  inside  the  Voronoi  cell  elements.  The  matrix  product  may  be  expressed  as 


[Glcr]  -[Gcr] 

< 

1 

q" 

2 

'  =  [Uj[Aj[V]  < 

1 

qcr 

2 

'  —  [U]  [A]  < 

1 

qcr 

2 

qcr 

k  4 

q^ 

.  4 

[G~]  -[dhll 


(4.73) 


[Gcrj  -[G" 


[A]  is  a  rectangular  matrix 


[U]  and  [V]  are  orthonormal  matrices  obtained  by  SVD  of 
with  nonnegative  values  on  the  diagonal.  The  zero  or  singular  (very  small  values  in  numerical  computations) 
values  in  [A]  corresponds  to  either  trivial  solutions  or  rigid  body  modes  of  the  displacement  solution.  For 
accurate  displacements,  elements  in  {q”-}  corresponding  to  small  or  zero  eigen-values  in  [A]  are  eliminated. 


4.6  Numerical  Examples 

The  numerical  examples  solved,  are  divided  into  four  categories.  In  the  first  set  of  examples,  the  convergence 
of  X-VCFEM  enriched  by  multi-resolution  wavelet  functions  is  demonstrated  for  static  cracks  by  comparison 
with  theoretical  predictions  and  results  available  in  the  literature.  The  second  set  of  examples  show  the 
effectiveness  of  X-VCFEM  in  modeling  the  propagation  of  multiple  cohesive  cracks.  The  third  set  of  examples 
is  intended  to  investigate  the  effect  of  cohesive  parameters  on  crack  growth.  The  final  set  of  examples  looks 
into  the  growth  of  multiple  pre-existing  cracks  to  comprehend  the  effect  of  morphology,  e.g.  distribution, 
orientation  etc.. 

4.6.1  Convergence  tests  for  X-VCFEM  for  static  cracks 

Effects  of  translation  and  dilation  parameters 

Figure  4.10(a)  shows  a  center  cracked  plate  of  width  2w=4cm  and  length  b=12cm  with  a  crack  length 
of  2a=  1.6cm.  The  plate  is  loaded  in  simple  tension  with  a  constant  remote  load  of  cr  =  5  MPa.  The 
material  parameters  are:  Young’s  modulus  E  =  1  MPa  and  Poisson  ratio  v  =  0.3.  Due  to  problem 
symmetry,  only  the  right  half  of  the  plate  is  modeled  with  one  X-VCFEM  element,  as  shown  in  figure 
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4.10(b).  Symmetry  conditions  are  imposed  on  the  left  edge.  The  crack  face  is  modeled  using  10  node-pairs 
and  the  element  boundary  consists  of  22  segments.  The  stress  function  used  in  this  example  consists  of 
the  three  parts  discussed  in  section  4.3.  For  the  polynomial  function,  the  order  of  interpolation  in  equation 
(4.19)  corresponds  to  pn  =  13  and  qn  =  13  for  a  total  of  102  terms.  For  the  branch  function  in  equation 
(4.23)  consists  of  only  1  term  with  sn  =  0  and  tn  =  0  .  The  wavelet  functions  are  changed  from  a  lower  level 
to  a  higher  level  using  the  adaptation  criterion  discussed  in  section  4.3.4.  Similar  parameters  are  assumed 
for  the  and  771  directions,  i.e.,  aj  —  a,  tra  =  trc  and  b  1  =  d,\  in  equation  (4.37).  The  starting  values  of 
the  parameters  for  the  lower  level  (m  =  k  =  1)  are:  —  2  <  n  <  3,  0  <  /  <  1,  and  ai  =  C\  =  0.1.  The 
result  of  X-VCFEM  for  this  problem  is  plotted  in  terms  of  the  stress  ayy  along  the  crack  face  ( y  =  0) 
as  a  function  of  the  distance  from  the  center  of  the  crack  in  figure  4.11.  The  figure  4.11(a)  corresponds  to 
the  stresses  by  varying  the  translation  parameter  Zq,  while  figure  4.11(b)  is  for  the  variation  of  the  dilation 
parameter  a\.  From  figure  4.11(a)  it  is  evident  that  a  smaller  b  1  make  the  stress  concentration  at  the  crack 
tip  higher.  However,  very  small  61  <  0.001  (no  translation)  leads  to  linear  dependence  of  the  columns  of 
the  [P]  matrix  generated  from  the  wavelet  basis  functions,  and  should  be  avoided.  Figure  4.11(b)  shows 
that  smaller  a\  results  in  faster  convergence  to  higher  crack  tip  stress  concentration.  However,  very  small 
values  of  a  1  can  also  lead  to  oscillatory  stresses.  On  the  other  hand,  large  01  values  («  0.15)  shifts  the 
stress  peak.  The  optimal  selection  of  these  parameters  is  therefore  very  important.  This  is  obtained  through 
the  multi-resolution  construction  of  bases,  discuss  next.  The  multi-resolution  wavelet  bases  are  significantly 
more  effective  in  simulating  crack  problems.  Table  4.7  shows  the  effect  of  the  dilation  transfer  rate  tra  =  trc 
on  the  stress  intensity  factors  for  variation  in  the  translation  parameters  61  =  d\.  Other  parameters  in  the 
simulation  are  a\  =  C\  =0.1,  mn  =  kn  =  3,  nn  =  6,  and  ln  —  2.  The  values  tra  =  trc  —  1  imply  no  dilation. 
As  tra  approaches  1,  the  different  levels  functions  become  more  and  more  dependent  on  each  other.  From 
the  table,  the  minimum  error  is  achieved  for  61  =0.1  and  tra  =  trc  =  0.5  or  0.6. 

Convergence  with  multi-resolution  wavelet  bases 

The  example  in  section  4.6.1  is  considered  again  for  studying  the  solution  convergence  behavior  with  multi¬ 
resolution  wavelet  functions.  The  four  sets  of  parameters  represent  four  instances  of  multi-resolution  stress 
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function  enrichment.  The  first  case  consists  of  only  polynomial  and  branch  functions  for  the  stress  inter¬ 
polation,  for  which  the  details  are  provided  in  section  4.6.1.  Cases  2,  3,  4,  and  5,  introduce  different  levels 
of  the  wavelet  basis  functions.  The  wavelet  parameters  common  to  these  four  cases  are:  nn  =  6 ,  /„  =  2, 
aj  =  Ci  =  &i  =  di  —  0.1,  and  tra  =  trc  =  0.5.  The  parameters  corresponding  to  the  levels  of  the  multi¬ 
resolution  enrichment  (mn  =  kn)  are  listed  in  table  4.7. 

The  mode  I  stress  intensity  factor  is  calculated  for  all  the  five  cases  and  is  normalized  with  respect  to  the 
analytical  prediction  Krej  by  linear  elastic  fracture  mechanics  (LEFM),  reported  in  [98].  The  second  row  of 
table  4.7  compares  this  value  for  the  different  cases.  Without  the  wavelet  bases  (case  1),  the  solution  is  16% 
higher  than  the  theoretical  value.  Cases  2-5  results  demonstrate  that  the  wavelet  basis  effectively  reduces 
the  error  with  increasing  resolution  level  (mn).  The  X-VCFEM  generated  stress  ayy  at  y  =  0  is  plotted  in 
figure  4.12  for  cases  1-4.  Without  the  wavelet  enrichment,  the  stress  concentration  at  the  crack  tip  ( x  =  0.8) 
is  completely  misrepresented.  The  stress  peaks  are  represented  with  increasing  accuracy  with  additional 
levels  of  multi-resolution  wavelet  functions.  The  strain  energy  error  in  equation  (4.45)  is  also  calculated  for 
the  cases  2-5  and  tabulated  in  table  4.7.  The  error  rapidly  decreases  with  increasing  wavelet  enrichment, 
confirming  the  fast  convergence  rate  of  the  multi-resolution  algorithm.  However,  the  stress  intensity  factor 
K[  is  calculated  from  a  contour  that  is  away  from  the  crack  tip.  The  stresses  on  this  contour  are  much 
more  stabilized  and  additional  wavelet  bases  do  not  affect  these  stresses  considerably.  Hence,  the  error  in 
Ki  is  not  significantly  affected  by  their  addition.  From  the  above  convergence  tests,  the  optimal  parameters 
for  stress  function  representations  in  X-VCFEM  are  chosen  to  be  pn  =  qn  =  13,  sn  =  tn  =  0,  nn  =  6, 
ln  =  1,  mn  =  kn  =  4,  a\  —  a  —  6i  =  d-i  =  0.1  and  tra  =  trc  =  0.5.  These  are  retained  for  all  subsequent 
simulations.  X-VCFEM  simulations  of  the  cracked  plate  are  further  conducted  for  different  crack  lengths, 
to  study  the  effect  of  this  length  on  the  solution  convergence.  The  specific  dimensions  in  figure  4.10(a) 
are  2w=2  cm,  b=6  cm,  while  the  crack  length  2a  is  varied.  The  plate  is  loaded  under  remote  tension  of 
a  —  40Pa.  X-VCFEM  solution  of  Kj  for  various  values  of  a/w  are  plotted  in  figure  4.13  and  compared 
with  the  theoretical  predictions  of  [98].  X-VCFEM  predictions  match  the  theoretical  results  extremely  well. 
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4.6.2  Efficiency  and  Accuracy  of  X-VCFEM 


Prior  to  studying  the  effect  of  cohesive  parameters  and  multi-crack  distributions,  the  accuracy  and  efficiency 
of  X-VCFE  model  are  validated  by  several  numerical  examples. 

Comparing  efficiency  with  ABAQUS  for  a  simple  crack  propagation  problem 

A  plate  with  a  pre-existing  edge  crack  under  remote  tension  load  is  solved  for  plane  strain  by  X-VCFEM 
and  ABAQUS  as  shown  in  figure  4.14(a).  The  material  Young’s  modulus  E  =  70,000 MPa,  and  Poisson 
ratio  v  =  0.33.  A  bilinear  cohesive  zone  model  discussed  in  [58]  is  used  to  describe  the  crack  growth  and 
the  cohesive  model  parameters  are  amo.x  =  5  M Pa,  6C  =  1  x  10~6  mm,  5e  =  5  x  10-3  mm  and  /?  =  0.707. 
The  entire  domain  is  represented  by  a  single  element  in  X-VCFEM,  consisting  of  142  nodes  for  displacement 
interpolation.  The  adaptive  enrichment  of  wavelet  bases  is  determined  by  the  strain  energy  error  in  equation 
(4.45).  As  shown  in  previous  section,  the  optimal  parameters  for  stress  function  representations  in  X-VCFEM 
are  chosen  to  be  pn  —  qn  —  13,  sn  =  tn  =  0,  nn  —  8,  ln  =  1,  mn  =  kn  —  4,  ai  =  Ci  =  b\  =  di  =  0.1  and 
tra  =  trc  =  0.5,  which  means  that  stress  function  interpolations  consist  of  102  terms  of  polynomial  functions, 
1  term  in  the  branch  function,  and  128  terms  in  the  wavelet  function  representation.  These  are  retained  for 
all  subsequent  simulations.  It  is  assumed  that  the  crack  propagates  horizontally  due  to  symmetry,  and  hence 
the  modules  for  determining  incremental  crack  direction  in  section  4.5  is  switched  off  for  this  problem.  A 
special  UEL  subroutine  is  developed  in  ABAQUS  for  incorporating  the  cohesive  model  at  a  given  interfaces. 
A  total  of  12840  4-node  2D  element  and  77  cohesive  elements  are  used  in  ABAQUS.  Figure  4.14(b)  shows 
the  load  cr-vertical  displacement  uy  plot  at  point  A.  The  two  codes  yield  very  similar  results,  an  attestation 
of  X-VCFEM  accuracy.  However,  the  X-VCFEM  simulation  takes  only  1.6  minutes  on  a  single  CPU  in  the 
Pentium  4  cluster  with  2.4Ghz  Intel  P4  Xeon  processors,  as  opposed  to  13.9  minutes  by  the  ABAQUS  run  on 
the  same  machine.  Thus,  even  for  this  simple  example,  a  tenfold  advantage  in  computing  speed  is  achieved 
by  X-VCFEM.  It  is  expected  that  this  factor  will  increase  considerably  with  increasing  complexity,  such  as 
more  cracks. 


85 


A  classical  problem  on  dynamic  crack  propagation 

This  numerical  example  is  based  on  Kalthoff’s  well  known  experiment  on  dynamic  crack  propagation  in  a 
impact  loaded  prenotched  plate,  that  has  been  the  subject  of  many  studies  [55,  56,  85].  These  studies  suggest 
that  a  crack,  subjected  to  a  tension-compression  load  as  shown  in  figure  4.15(a),  propagates  at  an  angle  of 
approximately  60°  —  70°  with  respect  to  the  initial  notch  in  the  plate.  The  present  X-VCFEM  does  not 
explicitly  incorporate  inertia  terms,  and  hence  a  quasi-static  crack  propagation  problem  is  simulated  instead 
of  the  dynamic  test.  The  configuration  in  figure  4.15(a),  shows  that  the  experimental  projectile  motion  is 
replaced  by  the  traction  boundary  conditions  in  the  simulation  under  plane  strain  conditions.  A  small  initial 
crack  length  of  o=0.02 m  is  chosen  to  mitigate  the  effect  of  the  constrained  right  hand  boundary  on  crack 
propagation.  Material  properties  for  this  problem  are:  Young’s  modulus  E  =  207  GPa  and  Poisson  ratio 

v  —  0.3  and  cohesive  zone  model  parameters  in  equation  (4.12)  are:  amax  =  0.1  MPa,  6e  —  lx  10~6m,  and 

P  =  0.  The  entire  domain  is  represented  by  a  single  element  in  X-VCFEM,  consisting  of  132  nodal  degrees 
of  freedom.  The  results  of  the  X-VCFEM  simulation  is  shown  in  figure  4.15.  From  figure  4.15(a)  the  initial 
crack  growth  angle  is  around  70°,  which  is  corroborated  by  brittle  failure  experiments  at  very  low  velocities 
[32],  Subsequently,  the  crack  propagation  takes  place  within  the  envelope  of  60°  —  70°,  which  is  in  agreement 
with  studies  in  [55,  56,  85].  The  dynamic  conditions,  as  well  as  boundary  constraints  are  responsible  for  the 
small  difference  between  X-VCFEM  results  and  those  in  [85].  Volume-averaged  or  macroscopic  shear  stress- 
shear  strain  behavior  for  this  problem  is  plotted  in  figure  4.15(b).  The  volume  averaging  of  the  local  stress 
and  strain  fields  over  the  entire  microscopic  domain  fi  is  performed  as 

d-.j(t)  =  ^  J^aij(xk,t)dV 

*ij(t)  =  ^  Jeij(xk,t)dV  (4.74) 

where  xk  and  t  are  the  spatial  coordinates  and  time  (cumulative  increments  in  these  problems)  respectively, 
and 

^  (Mt)K'  +  Mi)K)<idfi  (4.75) 
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aij(t )  represents  the  effective  strain  field  caused  by  the  possible  displacement  jump  at  the  crack.  It  is 
calculated  along  the  crack  path  rcr  with  [iit (i)]  denoting  the  displacement  jump. 

Crack  propagation  in  sheared  plate  with  a  central  crack 

This  example  is  based  on  a  classical  problem  of  a  single  crack  propagation  in  a  large  plate  with  a  central  crack. 
The  plate  is  subjected  to  a  far  field  shear  load.  The  problem  was  experimentally  studied  by  Erdogan  and  Sih 
[32]  and  an  optical  micrograph  of  their  cracked  specimen  is  shown  in  figure  4.16(a).  The  specimen  material  in 
their  experiment  were  assumed  to  be  homogeneous,  isotropic  and  linearly  elastic  and  the  crack  was  assumed 
to  be  brittle.  A  single  element  of  dimension  10  m  x  8  m  in  X-VCFEM  is  used  to  simulate  this  experiment  as 
shown  in  figure  4.16(b).  The  initial  crack  length  is  a=1.6  m.  The  material  parameters  are:  Young’s  modulus 
E  =  100  GPa,  Poisson  ratio  u  =  0.3  and  the  cohesive  law  parameters  are:  amax  =  0.1  MPa ,  (3  =  1,  and 
5C  =  1  x  10~7m.  The  shear  stress  applied  on  the  top  and  bottom  surfaces,  is  varied  from  0  to  0.041  GPa  with 
plane  stress  assumptions.  As  shown  in  figure  4.16(b),  the  crack  path  predicted  by  X-VCFEM  compares  well 
with  the  observations  in  [32],  Figure  4.16(c)  shows  the  growth  of  the  crack  opening  displacement  components 
at  the  right  tip  A.  The  entire  computational  process  took  20  minutes  on  a  single  CPU  in  the  Pentium  4 
cluster  with  2.4Ghz  intel  P4  Xeon  processors. 

Crack  propagation  in  three-point  bending  specimen 

Two  numerical  examples  are  considered  for  this  specimen.  In  the  first  example,  symmetric  mode  I  crack 
propagation  in  a  three-point  bending  test,  as  shown  in  figure  4.17,  is  modeled.  Plane  stress  conditions  are 
assumed  in  the  simulation.  This  problem  of  cohesive  crack  propagation  has  been  studied  by  Carpinteri  [17] 
using  node  release  technique  and  by  Moes  and  Belytschko  [66]  using  the  extended  FEM  or  XFEM.  The 
geometrical  dimensions  for  the  specimen  in  figure  4.17  are  b=  0.15  m,l=4b,  t(specimen  thickness)=  b,a=0, 
and  d=0.001  m.  The  material  properties  are:  Young’s  modulus  E  =  36,500  MPa,  Poisson  ratio  v  —  0.1, 
and  the  cohesive  parameters  are  amax  =  3.19  MPa,  and  P  —  0.  The  X-VCFEM  solution  is  compared 
with  that  in  [66]  through  the  load-deflection  curve  of  figure  4.18.  The  cohesive  displacement  parameters  are 
5e  =  3.134796  x  10~5  m  and  Se  =  6.26959  x  10~6m  for  figures  4.18(a)  and  4.18(b)  respectively.  A  sharper 
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snap-back  is  seen  for  the  latter  case.  Excellent  match  is  observed  between  the  X-VCFEM  and  XFEM  results 
The  second  example  shows  a  mixed-mode  cohesive  crack  propagation  in  a  three-point  bend  test  due  to  an 
unsymmetrically  positioned  initial  crack.  The  problem,  shown  in  figure  4.19(a),  has  been  studied  by  Mariani 
and  Perego[64]  using  XFEM  under  plane  stress  conditions.  The  initial  crack  position  is  determined  by  the 
offset  ratio  a,  defined  as  the  ratio  of  the  distance  of  the  initial  crack  from  the  mid-span  cross-section  to  half 
of  the  beam  span.  The  material  Young’s  modulus  E  =  31370  MPa,  and  Poisson  ratio  u  =  0.2.  The  cohesive 
model  parameters  are  <Jmax  =  4.4  MPa,  Se  =  0.07719298  mm  and  p  =  1.0.  Once  again,  the  entire  domain 
is  represented  by  a  single  element  in  X-VCEFM  with  154  nodal  degrees  of  freedom.  Figure  4.19(b)  shows 
the  load-deflection  curve  for  two  values  of  the  offset  parameter,  i.e.  a  —  0.5  and  a  =  0.25.  The  initial  elastic 
response  in  the  load  P-displacement  u  curve  is  stiffer  and  also  the  peak  load  is  higher  for  higher  values  of 
a.  The  load-displacement  response  exhibits  softening  in  the  later  stages  of  crack  propagation  due  to  the 
significantly  evolved  crack.  The  path  of  crack  propagation  for  the  two  cases  are  shown  in  figures  4.19(c)  and 
(d).  The  cracks  move  towards  the  point  of  applied  load  and  align  themselves  perpendicular  to  the  edge  of 
the  specimen. Excellent  agreement  is  obtained  between  the  results  by  X-VCFEM  and  in  [64], 

4.6.3  Mesh  independence  of  crack  propagation  with  X-VCFEM 

A  panel  with  domain  5  cm  x  3  cm  containing  two  initial  cracks  is  remotely  loaded  in  tension  as  shown  in 
figure  4.20(a).  The  problem  has  been  solved  by  Sharma  et.  al.  [95]  using  the  element  free  Galerkin  meshless 
method.  For  X-VCFEM  solution,  the  domain  is  meshed  into  two  elements  with  three  different  topologies 
shown  in  figure  4.20.  Plane  stress  conditions  are  again  assumed.  A  total  of  11  increments  is  used  to  model 
the  entire  crack  propagation  process.  The  material  parameters  are:  Young’s  modulus  E  =  207  GPa  and 
Poisson  ratio  u  =  0.3  and  cohesive  zone  parameters  are:  amax  =0.1  MPa,  6e  =  1  x  10-6cm,  and  p  =  1. 
The  three  figures  4,20(b,c,d)  show  no  mesh  dependence  of  the  X-VCFEM  predictions  and  the  comparison 
with  results  in  [95]  is  excellent. 
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4.6.4  Effect  of  cohesive  parameters  on  crack  evolution 


This  example  is  intended  to  investigate  the  effect  of  cohesive  parameters  on  crack  growth.  Cohesive  zone 
model  parameters,  e.g.  crmax  and  5e  in  equation  (4.12),  can  significantly  affect  the  propagation  and  overall 
behavior  of  a  cracking  material.  The  effects  of  these  cohesive  parameters  are  studied  for  crack  propagation 
in  a  sheared  plate  with  a  central  crack  subjected  to  a  far  field  shear  load.  This  classical  problem  was 
experimentally  studied  by  Erdogan  and  Sih  [32]  and  an  optical  micrograph  of  their  cracked  specimen  is  shown 
in  figure  4.21(a).  The  specimen  material  in  their  experiment  was  assumed  to  be  homogeneous,  isotropic  and 
linearly  elastic  and  the  crack  was  assumed  to  be  brittle.  A  single  element  of  dimension  10  m  x  8  m  in  X- 
VCFEM  is  used  to  simulate  this  experiment  as  shown  in  figure  4.21(b).  The  initial  crack  length  is  lo=l-6  m. 
The  material  parameters  are:  Young’s  modulus  E  =  100  GPa ,  Poisson  ratio  v  =  0.3.  Five  different  sets  of 
cohesive  parameters,  illustrated  in  figure  4.21(b)  are  considered  for  this  example.  These  are 

•  A:  o-max=3.0  MPa,  <5e =3.0  e-4m,  p  =  1.0 

•  B:  (jmax=6.0  MPa,  <5e  =  1 . 5  e-4m,  p  =  1.0 

•  C:  crmQX=3.0  MPa,  <5e  =6.0  e-4m,  p  =  1.0 

•  D:  amax—  6-0  MPa,  <5e=3.0  e-4m,  p  =  1.0 

•  E:  <rmax  =  1.5  MPa,  <5e=6.0  e-4m,  p  =  1.0 

As  shown  in  figure  4.21(b),  all  the  cases  correspond  to  the  same  cohesive  energy.  The  load  is  applied  by 
controlling  the  opening  of  crack  propagation  through  fixed  values  of  the  increment  A l  in  equation  (4.55). 
Further  a  uniform  shear  load  per  unit  length  to  is  applied  on  the  top  and  bottom  surfaces  as  shown  in  figure 
4.21(c).  In  each  increment,  the  applied  load  is  scaled  by  the  arc-length  parameter  A  of  equation  (4.54),  to 
yield  an  equilibriated  applied  load  corresponding  to  a  prescribed  crack  propagation  length.  The  crack  path 
for  all  the  different  cohesive  parameters  predicted  by  X-VCFEM  are  very  similar  and  compare  well  with 
experimental  observations  in  [32],  However  a  considerable  dependence  on  cohesive  parameters  is  seen  in  the 
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shear-crack  length  response,  demonstrated  in  figure  4.21(d),  where  the  normalized  crack  length  is  defined  as 


The  current  crack  length 
The  initial  crack  length 


(4.76) 


This  points  to  the  fact  that  the  rate  of  propagation,  and  not  the  direction,  is  dependent  on  the  parameters 
for  this  problem.  For  the  cases  with  larger  peak  traction  cases:  B  and  D,  higher  applied  loads  are  needed 
for  causing  similar  crack  growths  as  for  cases  with  lower  peak  traction:  A  and  C.  Comparison  of  the  results 
for  cases  B  and  D,  show  that  a  smaller  5e  (case  B)  results  in  quicker  reduction  of  the  local  cohesive  traction. 
This  makes  the  overall  load  for  the  case  B  to  increase  slower  than  that  for  case  D  with  a  higher  8e.  The  case 
E  consistent  with  the  trends  exhibited  by  the  other  load  cases.  Although,  the  simulation  results  show  that 
both  amax  and  5e  affect  the  crack  growth,  comparison  of  cases  A,  B,  C  with  D  shows  that  the  crack  growth 
is  more  sensitive  to  (Tmax  than  to  <5e.  The  results  also  imply  that  the  cohesive  energy,  or  effectively  the  energy 
release  rate  Gc  does  not  alone  determine  the  properties  of  crack  propagation.  The  individual  parameters, 
affecting  the  shape  of  the  cohesive  law,  play  an  important  role  in  predicting  the  growth  characteristics.  These 
effects  are  also  tested  for  multiple  crack  growth  in  the  next  set  of  examples. 


4.6.5  Propagation  of  multiple  pre-existing  cracks 

The  final  set  of  examples  looks  into  the  growth  of  multiple  pre-existing  cracks  to  comprehend  the  effect  of 
morphology,  e.g.  distribution,  orientation  etc.. 

Firstly,  a  plate  with  five  randomly  located  cracks  is  simulated  under  a  tensile  loading  as  shown  in  figure 
4.22(a).  The  plate  has  dimensions  0.6  m  x  0.4  m;  material  parameters:  Young’s  modulus  E  =  105  MPa 
and  Poisson  ratio  v  =  0.3;  and  cohesive  parameters:  amax  =  0.1  MPa,  0=1,  and  5e  =  lx  10~5cm.  Figure 
4.22(b)  shows  the  final  positions  of  the  cracks  that  have  grown  with  the  loading.  The  cracks  propagate  across 
element  boundaries  and  are  attracted  to  each  other  in  certain  regions  till  they  nearly  merge. 

A  plate  with  28  randomly  located  and  oriented  cracks  is  simulated  under  a  tensile  loading.  Figures  4.23(a) 
and  (b)  show  the  two  microstructures  with  different  crack  distributions.  For  the  microstructure  1,  all  the 
cracks  of  equal  length  are  oriented  horizontally  and  their  distribution  is  random.  The  microstructure  2  has 
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cracks  of  random  length  and  orientation.  In  addition,  it  contains  a  cluster  of  8  cracks  in  a  otherwise  random 
distribution  as  shown  in  figure  4.23(b).  The  plate  is  of  dimension  O.lmxO.l  m,  and  the  material  parameters 
are:  Young’s  modulus  E  =  104  MPa  and  Poisson  ratio  u  =  0.3.  To  understand  the  effect  of  cohesive 
parameters  on  crack  propagation,  two  different  sets  of  cohesive  parameters  are  considered.  They  are: 

CP-1:  crmax=l.O  MPa,  $e  =  1.0  e-5m,  (3  —  0.707 
CP-2:  amax=2.0  MPa,  Se= 0.5  e-5m,  0  =  0.707 

A  uniform  tension  load  per  unit  length  a  is  applied  on  the  top  and  bottom  surfaces  as  shown  in  figure 
4.23(a,b).  In  each  increment,  the  applied  load  is  scaled  by  the  arc-length  parameter  A  of  equation  (4.54),  to 
yield  an  equilibriated  applied  load  corresponding  to  a  prescribed  crack  opening  deformation. 

Figures  4.24(a,b)  and  (c,d)  show  the  contour  plots  of  the  microstructural  stress  ayy  together  with  evolved 
position  of  the  cracks  at  the  final  stage  of  loading,  for  the  two  sets  of  microstructures  and  cohesive  parameter 
respectively.  The  growth  pattern  of  each  crack  can  be  observed  by  comparing  with  its  initial  configuration 
in  figures  4.23(a)  and  (b).  The  cracks  propagate  across  element  boundaries,  interact  with  each  other  and  in 
some  cases,  they  merge.  The  relation  of  the  propagation  of  multiple  cracks  to  the  morphology  and  cohesive 
parameters  is  in  general  complicated.  However,  several  observations  can  be  made  based  on  the  results  of  the 
simulation  by  this  model. 

•  Larger  stress  concentrations  develop  at  tips  of  cracks  that  are  nearly  perpendicular  to  the  direction  of 
loading.  Consequently,  this  subset  of  cracks  grow  more  easily  than  others  that  are  more  aligned  with 
the  loading  direction.  From  figure  4.24(b)  and  (d),  it  can  be  seen  that  some  cracks  that  are  nearly 
parallel  to  the  load  direction  never  propagate. 

•  Stress  concentrations  are  higher  at  tips  of  longer  cracks.  The  reason  is  stress  concentrations  at  crack 
tips  come  from  the  external  load,  which  cannot  be  handled  by  the  weak  crack.  Longer  cracks  lead  to 
more  external  load  concentrating  to  tips.  This  is  verified  by  results  shown  in  figure  4.24(b)  and  (d), 
where  longer  cracks  are  easier  to  propagate  than  shorter  ones. 

•  Irrespective  of  the  initial  orientation,  the  evolved  crack  path  tends  to  align  in  a  direction  perpendicular 
to  the  applied  load  direction.  This  correspond  to  an  optimal  direction  for  releasing  the  cohesive  energy. 
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This  observation  is  dominant,  when  the  influence  of  nearby  cracks  on  the  local  stress  field  is  small. 
The  local  stress  field  for  this  phenomenon  is  mainly  governed  by  the  influence  of  applied  load  on  this 
single  crack. 

•  Cracks  are  attracted  towards  weak  surfaces,  such  as  other  cracks  or  voids  and  prefer  to  propagate  in 
those  directions.  This  may  be  attributed  to  the  fact  that  the  cohesive  energy  in  the  direction  of  these 
weaker  surfaces  with  lower  (or  zero)  tractions  is  naturally  lower  in  comparison  with  other  directions. 

•  Figures  4.24(b)  and  (d)  show  that  the  longest  crack  does  not  necessarily  evolve  from  a  cluster.  Not  all 
cracks  in  a  cluster  grow  considerably.  This  is  somewhat  in  contrast  to  observations  made  with  particle 
reinforced  composites,  where  almost  always  clusters  cause  a  local  stress  concentration.  The  interaction 
between  neighboring  cracks  contributes  to  the  enhancement  or  mitigation  of  stresses,  depending  on 
their  orientations  and  length.  This  dictates  their  propagation,  and  just  being  in  a  cluster  does  not 
guarantee  significant  growth. 

•  The  different  cohesive  parameters  show  very  little  difference  in  the  final  configuration  and  hence  the 
propagation  direction.  However,  the  rate  of  crack  growth  varies  considerably  with  these  parameters  as 
seen  in  the  crack  length-macroscopic  strain  plot  of  figure  4.25 

Figure  4.26  shows  the  macroscopic  stress-strain  response  for  the  two  microstructures  and  cohesive  pa¬ 
rameters.  Even  before  the  cracks  propagate  (corresponding  to  the  change  in  slope),  the  stiffness  of  the 
micro  structure  2  is  higher  than  that  of  micro  structure  1  due  to  a  higher  level  of  effective  damage  caused  by 
crack  lengths  and  more  importantly  orientations.  Orientations  perpendicular  to  the  load  direction  causes  a 
larger  reduction  in  stiffness  in  comparison  with  other  directions.  With  additional  loading,  the  overall  damage 
caused  by  the  growth  of  cracks  is  also  higher  for  the  microstructure  1.  This  is  seen  by  the  lower  values  of  the 
macroscopic  stress  for  this  case.  The  effect  of  the  cohesive  parameters  on  the  stress-strain  response  is  quite 
pronounced.  The  maximum  macroscopic  stress  for  both  microstructures  increases  significantly  for  higher 
values  of  <7max,  even  though  the  cohesive  energy  is  the  same  for  the  two  cohesive  models.  This  is  caused  by 
a  slowdown  in  the  growth  rate  of  the  cracks  with  overall  deformation. 
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4.7  Concluding  Remarks 


The  extended  Voronoi  cell  finite  element  model  is  developed  in  this  chapter  to  predict  initiation  and  growth 
of  damage  by  crack  propagation  in  brittle  matrix.  The  cracks  are  modeled  by  a  linear  cohesive  zone  model 
and  their  incremental  directions  and  growth  lengths  are  determined  in  terms  of  the  cohesive  energy  near  the 
crack  tip.  Important  enhancements  are  made  to  the  element  to  allow  stress  discontinuities  across  the  cohesive 
crack  and  to  accurately  depict  the  crack  tip  stress  concentrations.  These  features  are  accommodated  through 
the  incorporation  of  (a)  branch  functions  in  conjunction  with  level  set  methods  across  crack  contours,  and 
(b)  adaptive  multi-resolution  wavelet  functions  in  the  vicinity  of  the  crack  tip.  Several  problems  are  solved 
and  compared  with  existing  solutions  in  the  literature  for  validation  of  the  X-VCFEM  algorithms,  both  with 
respect  to  macroscopic  (load-deformation  behavior)  and  microscopic  (crack  path).  The  X-VCFEM  results 
show  excellent  accuracy  in  their  comparison  with  analytical  and  other  numerical  solutions.  Also  compari¬ 
son  with  ABAQUS  shows  the  efficiency  of  X-VCFEM.  Numerical  simulations  are  conducted  with  different 
<7max  and  <5C  to  understand  the  effect  of  cohesive  parameters  on  the  crack  propagation.  It’s  observed  that  in 
addition  to  the  total  cohesive  energy,  the  individual  parameters  have  effects  on  crack  growth.  The  effect  of 
geometrical  information  of  multiple  pre-existing  cracks,  including  the  lengths,  positions  and  orientations  of 
cracks,  on  their  propagation  is  studied  by  simulating  a  plate  with  28  randomly  located  and  oriented  cracks. 
Simulation  results  show  that  the  crack  with  a  longer  length  and  nearly  perpendicular  to  load  direction  is 
easier  to  propagation  than  other  cracks.  Cracks  propagation  direction  is  dependent  on  the  local  stress  field, 
which  is  managed  by  both  the  external  load  and  nearby  material  phases,  such  as  other  cracks  in  a  cluster. 
This  research  reveals  the  significance  of  analyzing  large  regions  of  the  microstructure  and  proves  the  effec¬ 
tiveness  of  the  X-VCFEM.  The  simulation  results  based  on  X-VCFEM  could  also  provide  positive  feedback 
for  design  modification. 

Based  on  the  study  on  interfacial  debonding  and  cohesive  matrix  cracking  in  composites,  the  interaction 
of  the  two  damage  phenomena  is  studied  in  the  next  chapter,  where  the  X-VCFEM  is  improved  and  a 
criterion  for  assessing  the  direction  of  damage  development  is  proposed. 
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tra 

0.6 

0.6 

0.6 

0.5 

0.5 

0.5 

0.4 

0.4 

0.4 

bi 

0.05 

0.1 

0.2 

0.05 

0.1 

0.2 

0.05 

0.1 

0.2 

Ki/Kref 

1.019 

1.014 

1.040 

1.017 

1.014 

1.027 

1.020 

1.020 

1.035 

Table  4.1:  Normalized  stress  intensity  factors  ( Ki/Kref )  for  different  values  of  tra  and  6]  in  the  multi¬ 
resolution  wavelet  representation. 


Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

TYln  —  kn 

0 

1 

2 

3 

4 

K,/KreI 

1.1642 

1.0361 

1.0208 

1.0062 

1.0020 

ASE 

96.45% 

45.91% 

7.06% 

3.01% 

Table  4.2:  Errors  with  varying  enrichment  order  of  multi-resolution  wavelet  functions  for  the  different  cases. 
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Figure  4.3:  (a)  A  schematic  diagram  of  a  crack  surface  showing  parameters  related  to  the  distance  functions; 
(b)  depiction  of  the  branched  stress  function  <f,!,ranc/l  near  a  crack  for  s  =  0,  t  =  0. 
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VCFEM 


Regions  covered  by  wavelet 
basis  functions 


Figure  4.4:  Distribution  of  multi-resolution  wavelet  bases  around  a  crack  tip:  (a)  Crosses  (x)  refer  to  the 
location  of  the  origin  of  the  basis  vectors  at  a  lower  level  corresponding  to  dilation  parameters  ( tra  and  trc) 
and  (b)  adaptively  upgraded  to  higher  level  wavelet  bases  with  the  addition  of  the  next  level  of  bases  at 
locations  indicated  by  the  (□). 
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Figure  4.5:  Load-deflection  behavior  in  a  3-point  bend  test  with  a  crack,  showing  the  softening  snap  back 
phenomenon. 
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Figure  4.7:  The  stress  axx  at  x  =  —0.3  for  a  double  cantilever  beam  to  demonstrate  the  effect  of  the 
branched  stress  function. 


(a)  (b) 


(<0 


Figure  4.8:  Algorithms  for  incremental  propagation  of  cohesive  cracks:  (a)  for  direction  and  incremental 
length,  (b)  a  cohesive  crack  going  through  the  inter-element  boundary,  (c)  for  merger  with  other  cracks  and 
(d)  for  evaluation  of  J—  integral. 
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Figure  4.9:  Subdivision  of  the  Voronoi  cell  element  for  Gaussian  quadrature,  with  a  higher  density  of 
integration  points  near  the  crack  tip. 
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(a)  (b) 


Figure  4.10:  (a)  A  center  cracked  plate  loaded  in  tension,  (b)  a  single  X-VCFEM  element  with  prescribed 
boundary  conditions 
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uy  at  point  A  (m) 


(b) 

Figure  4.14:  (a)  A  plate  with  an  edge  crack  under  remote  tension  load,  (b)  comparison  of  load-deformation 
curves  by  X-VCFEM  and  ABAQUS. 
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(b) 


Figure  4.15:  (a)  Prediction  of  the  crack  path  by  X-VCFEM  for  the  Kalthoff  experiment,  (b)  the  macroscopic 
stress-strain  response. 
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(c) 

Figure  4.16:  (a)  Optical  micrograph  showing  the  path  of  cracking  in  a  plate  with  a  central  crack  sub¬ 
jected  to  far-field  shear  [32],  (b)  corresponding  crack  crack  path  generated  by  X-VCFEM,  (c)  crack  opening 
displacement  at  the  tip  A. 


109 


Figure  4.17:  A  three-point  symmetric  bending  specimen. 
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(a) 


(b) 


Figure  4.18:  Comparison  of  normalized  load-deflection  curves  for  the  three-point  bending  beam:  (a)  Se  = 
3.134796  x  10~5  m  and  (b)  6e  =  6.26959  x  1CT6  m. 


Ill 


Figure  4.19:  (a)  A  three-point  bending  specimen  with  an  unsymmetric  initial  crack,  (b)  comparison  of 
load-deflection  curves  from  X-VCFEM  and  literature  [64]  ,  (c)  and  (d)  comparison  of  the  crack  paths  by 
X-VCFEM  with  that  in  [64]  for  a  =  0.25  and  a  =  0.5,  respectively. 


(b) 


Figure  4.20:  A  plate  with  two  cracks  in  arbitrary  locations  modeled  by  X-VCFEM  using  elements  of  different 
topologies  located  cracks,  (b,c  and  d)  show  crack  path  at  the  end  of  the  loading  for  the  different  elements 
and  also  a  comparison  with  [95]. 
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(c) 


(d) 


Figure  4.21:  (a)  Optical  micrograph  showing  the  path  of  cracking  in  a  plate  with  a  central  crack  subjected  to 
far-field  shear  [32],  (b)  5  different  sets  of  cohesive  parameters  for  X-VCFEM  simulations,  (c)  corresponding 
crack  path  generated  by  X-VCFEM,  (d)  comparison  of  the  growth  of  cracks  for  different  cases. 


Figure  4.22:  (a)  X-VCFEM  mesh  for  a  plate  with  five  randomly  located  cracks,  (b)  crack  paths  at  the  end 
of  loading 


Figure  4.23:  Crack  propagation  in  two  square  regions  containing  28  cracks  by  X-VCFEM:  (a)  domain  with 
horizontal  cracks  of  equal  length  and  random  distribution,  (b)  domain  with  random  orientation,  length  and 
distribution  of  cracks  but  containing  a  cluster. 
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Figure  4.24:  Crack  propagation  in  two  square  regions  containing  28  cracks  by  X-VCFEM:  (a,b)  contour  plots 
of  a yy  (MPa)  with  cohesive  parameters  CP-1  for  the  domains  in  figure  23  (a)  and  (b),  (c,d)  contour  plots  of 
cryy  (MPa)  with  cohesive  parameters  CP-2  for  the  domains  in  figure  23  (a)  and  (b) 
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Figure  4.25:  Comparison  of  the  growth  of  crack  A  in  micro  structure  1  with  different  cohesive  parameters. 
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rameters. 
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Chapter  5 


Extended  Voronoi  Cell  Finite  Element 
for  Modeling  Interfacial  Debonding 
with  Matrix  Cohesive  Cracking  in 
Fiber  Reinforced  Composites 

5.1  Introduction 

Interfacial  debonding  and  cohesive  cracks  propagation  in  brittle  matrix  are  two  important  damage  phenomena 
in  fiber-matrix  composites.  Experiments  show  that  the  two  damage  phenomena  appear  in  the  same  material, 
where  the  failure  often  starts  from  the  interface  between  fiber  and  matrix,  and  is  subsequently  advanced  into 
matrix.  Researches  regarding  a  crack  meeting  a  bimaterial  interface  to  either  deflect  along  the  interface  or 
penetrate  into  the  next  layer  were  made  in  [3,  45,  46,  63],  where  the  criterion  of  deflection  versus  penetration 
was  established  based  on  the  energy  release  rate  and  fracture  energy.  However,  the  present  research  is  aimed 
at  only  elastic  cases,  which  requires  that  the  fracture  process  zone  at  the  crack  tip  is  small  compared  to  the 
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size  of  the  crack  and  the  size  of  the  specimen.  In  chapters  3  and  4,  cohesive  zone  models  are  introduced  into 
VCFEM  to  study  damage  of  interface  and  matrix,  where  the  stress  field  in  composites  are  described  by  a 
set  of  specifical  functions  accurately.  And  the  effect  of  cohesive  parameters  and  morphological  distributions 
are  studied  as  important  factors  to  the  damage  process.  All  work  in  the  two  chapters  are  theoretical  and 
computational  preparation  for  solving  interfacial  debonding  problems  coupled  with  matrix  cracking. 

In  this  chapter,  a  criterion  based  on  cohesive  zone  models  is  proposed  for  assessing  the  direction  of  damage 
development.  The  improved  X-VCFEM  is  developed  for  modeling  both  the  growth  of  interfacial  debonding 
and  the  propagation  of  multiple  cohesive  cracks  in  the  brittle  matrix  of  fiber-reinforced  composites.  The 
mechanics  theories  and  numerical  algorithm  in  previous  chapters  are  organized  as  an  organic  whole,  not  just 
a  simple  superposition.  It  begins  with  the  X-VCFEM  formulation  and  numerical  implementation,  followed 
by  the  numerical  example  showing  the  effectiveness  of  this  model  and  the  interaction  of  interface  and  crack 
propagation. 


5.2  Extended  Voronoi  cell  FEM  formulation  for  composites  with 
interfacial  debonding  and  matrix  cracking 

The  Voronoi  cell  finite  element  mesh  for  a  microstructure  with  both  debonded  interfaces  and  cohesive  cracks 
is  shown  in  figure  5.1(a),  where  the  region  is  divided  into  an  unstructured  finite  element  mesh  of  arbitrary 
Voronoi  cells.  A  typical  Voronoi  cell  element  Ue  is  shown  in  figure  5.1  (b).  Each  VC  element  is  composed 
of  the  matrix  phase  the  inclusion  phase(ftc),  the  interface  (f 2;n),  and  cracks  (ficr),  such  that  fle  = 

fim  U^c  U  U  ^cr,  where  interface  and  cracks  are  consider  as  zero  thickness  regions.  The  element  outer 
boundary  consists  of  the  prescribed  displacement  boundary  (rum),  prescribed  traction  boundary  (rtm)  and 
the  inter-element  boundary  (rm),  so  i.e.  =  Tum  u  rtm  Urm.  Compatible  displacement  conditions  apply 
on  dflc.  dflrc  has  an  outward  normal  nc  (=nm),  while  ne  is  the  outward  normal  to  dfle.  In  order  to  describe 
debonding  with  progressing  deformation  through  decohesion,  the  interface  is  lined  with  a  set  of  node-pairs 
with  nodes  belonging  to  the  matrix  interface  (dfl™)  and  inclusion  interface  (dQ.cc)  respectively.  The  traction 
tcoh  between  node-pairs  on  the  crack  surface  are  modeled  by  the  cohesive  zone  traction-separation  law.  The 
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behaviors  of  cohesive  cracks  in  the  brittle  matrix  are  described  by  a  similar  method,  where  nodes  in  node-pairs 

1  2 

are  arranged  at  different  sides  of  a  crack  (rcr  and  rcr).  In  the  incremental  assumed  stress  hybrid  X-VCFEM 
formulation,  the  complementary  energy  functional  for  each  element  is  expressed  in  terms  of  increments  of 
stress  and  displacement  fields  as: 


ne (<x,  Act,  u,  Au) 
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(5.1) 


Here  B  is  the  complementary  energy  density  and  the  superscripts  m  and  c  correspond  to  variables  asso¬ 


ciated  with  the  matrix  and  inclusion  phases.  crm  and  erc  are  the  equilibrated  stress  fields,  em  and  ec  the 

1  2 

corresponding  strain  fields  in  different  phases  of  each  Voronoi  element.  Also,  ue,  um,  uc,  ucr  and  ucr  are 

1  2 

the  kinematically  admissible  displacement  fields  on  dfle,  dfl™,  dfl^,  rcr  and  rcr  respectively.  The  prefix  A 

corresponds  to  increments.  The  term  in  the  box  in  equation  (5.1)  provide  the  work  done  by  the  interfacial 

tractions  Tm  =  T™ nm  +  Ttm tm  due  to  interfacial  separation  (um  —  uc),  where  T™  and  Ttm  are  the  normal 

and  tangential  components  that  are  described  by  cohesive  laws  at  the  interface.  Similarly,  the  last  term 

provide  the  work  done  by  the  cohesive  tractions  T”-  =  T"ncr  +  T[rtcr  due  to  displacement  separation 
1  2 

(ucr  -  uCT)  along  the  crack,  where  T^r  and  Ttcr  are  the  normal  and  tangential  components  of  the  cohesive 
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force.  The  total  energy  for  the  entire  composite  domain  is  obtained  by  adding  the  energy  functionals  for  N 
elements  as 

N 

n  =  £ne  (5.2) 

e=l 

5.2.1  General  element  assumptions  and  weak  form 

In  the  absence  of  body  forces,  two  dimensional  stress  fields  satisfying  equilibrium  relations  can  be  generated 
from  the  Airy’s  stress  function  <F( x,y ).  In  the  incremental  formulation,  stress  increments  in  matrix  and 
inclusion  are  obtained  from  derivatives  of  the  stress  functions  A$m(a ;,  y)  and  A<J>c(a;,y)  as: 
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(5.3) 


where  {A/?m}  and  {A/3C}  are  the  column  of  unknown  stress  increment  coefficients.  Convergence  properties 
and  efficiency  of  X-VCFEM  depend  on  the  choice  of  <FTn.  These  functions  should  adequately  account  for  the 
geometry  and  location  of  the  heterogeneity  in  the  element,  so  stress  functions  for  matrix  are  decomposed 
into  (a)  a  purely  polynomial  function  (b)  a  reciprocal  function  $>™c,  (c)  a  branch  function  &franch 

and  (d)  wavelet  functions  $™vit  =  |v  +  $mc  + 

^branch  +  &wvu)'  The  selection  of  stress  functions  are 
discussed  in  chapters  3  and  4  detailed.  Inclusion  stress  functions  are  admitted  as  polynomial  function  <I>£oiy 
(*c  =  Compatible  displacement  fields  satisfying  inter-element  continuity  on  the  element  boundary 

dflf  and  intra-element  continuity  on  both  the  interface  and  the  crack  face  Fcr  are  generated  by 
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interpolation  of  nodal  displacements  as: 


{Aue}  =  [Le]{Aqe}  ondQe 
{Aum}  =  [Lm]{Aqm}  onSfi™ 

{Auc}  =  [Lc]{Aqc}  on  d£fc 
{A^}  =  [Lcr]{Aicr}  on  r‘r, 

{Aucr}  =  [Lcr){Aqcr}  on  rcr  (5.4) 

1  2 

The  interpolation  matrices  [Le],  [Lm],  [Lc],  [Lcr],  [ticr]  for  the  nodal  displacements  on  the  respective  bound¬ 
aries  are  constructed  using  standard  linear  or  hierarchical  shape  functions.  Since  nodes  on  the  inter¬ 
face  and  crack  surfaces  are  always  belonging  to  some  node-pair,  the  interpolation  matrices  are  chosen  as 

[Lm]  =  [I/]  and  [L"]  =  [Lcr], 

Substituting  the  interpolations  of  stress  and  displacement  fields  from  equations  (5.3)  and  (5.4)  into  equation 
(5.1)  results  in  the  matrix  form  of  the  element  complimentary  energy 


ne  =  -^{Pm  +  Apm}T{Hrn}{pm+Apm}-~{pc  +  Apc}T[Hc}{f3c  +  A/3c} 
+  {(3m  +  A^J^GHq6  +  Aqe}  -  {Pm  +  A/3m}7’[G]m{qm  +  Aqm} 

+  {(3C  +  A/3c}T[G]c{qc  +  Aqc}  +  {/3m  +  A/3m}T[Gcr]{qlcr  +  Aqcr} 


{(3m  +  A/3m}'  [G^Kq0"  +  Aqcr}  -  {t}r{qe  +  Aqc} 

p  ,(um+Aum-uc-Auc) 

-  /  Tm  •  d{ um  -  u c)ddQ 

Jd np/anj  um-ue) 

r  i'u*r+AuCP-u",-Au*r  !  2 

-  /  /,  3  Tcr  •  d(ucr  —  ucr)drcT. 

Jrcr  Jucr-ucr 


(5.5) 
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where 


[Hm]  =  [  [Pm]r(Sm][Pm]dftm  [Hc]  =  /  [Pc]T[Sc][Pc]dftc 
JQ'  J  ne 

[Ge]  =  /  [Pm]T[ne][Le](ian7n  [Gm]  =  /  [Pm]T[nm][Lm]dc>ftm 

Janc 

[Gc]  =  /  [Pc]T(nc][Lc]d30c  [GCT]  =  A  [P]T[ncr][Lcr]rfrcr 
Janc  J  rcr 

[Gcr]  =  /  [P]T[ncr][Lcr]rfrcr  {t}  =  /  [Le]T{t  + At}drtm 

3rv  drtm 


(5.6) 


Construction  of  appropriate  stress  functions  with  optimally  high  resolution  is  necessary  for  accurately  de¬ 
picting  high  stress  gradients  near  the  crack  tip. 


5.2.2  Solution  Method 


Crack  growth  in  multiply  cracked  materials  is  solved  using  an  incremental  approach,  where  a  set  of  elemental 
and  global  equations  are  solved  in  each  increment  for  stresses  and  displacements. 

1.  Local  equations  for  each  element  are  obtained  by  setting  the  variation  of  equation  (5.5)  with  respect 
to  the  stress  coefficients  A /3m  and  A/3C  to  zero.  This  results  in  the  weak  form  of  the  element  kinematic 
relations 
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1  2 

[Ge]  -[Gm]  [0]  [Gcr]  -[Gcr] 

[0]  [0]  [Gc]  [0]  [0] 


qm  +  Aqm 
qc  +  Aqc  ’ 

qlcr  +  Aqlcr 

2  2 
qcr  +  Aqcr 


or  in  a  condensed  form 


[H]e{/3  +  A/3}  e  —  [G]e{q  +  Aq}e 


(5.7) 


(5.8) 
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Since  equation  (5.8)  is  linear,  the  stress  coefficients  can  be  directly  expressed  in  terms  of  the  nodal  displace¬ 
ments,  provided  the  element  [Hje  matrix  is  invertible. 

2.  Subsequently,  the  weak  forms  of  the  global  traction  continuity  conditions  are  solved  by  setting  the  varia¬ 
tion  of  the  total  domain  energy  functional  with  respect  to  the  generalized  displacement  components  to  zero. 
This  results  in  the  weak  form  of  the  traction  reciprocity  conditions 
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(5.9) 


or  in  a  condensed  form: 


N 


y»r  = E^-je 

e=l  e=l 

The  forces  at  the  interface  and  crack  surface  are  expressed  in  terms  based  on  the  cohesive  energy  as 


(5.10) 
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Combining  equations  (5.8)  and  (5.10)  and  eliminating  the  stress  coefficients  {/3+A/3}e,  results  in  the  equation 
for  solving  the  generalized  displacements 


E{[G]nH]-[G]e}{q  +  Aq}  =  £{Telt}e  (5.13) 

e=l  e=l 

Equation  (5.13)  is  a  nonlinear  matrix  equation  system  due  to  the  cohesive  laws.  Consequently,  a  Newton- 
Raphson  iterative  solver  is  invoked  to  solve  for  the  increments  of  nodal  displacements.  The  linearized  form 
of  equation  (5.13)  for  the  y-th  iteration  is 

( £{Teit}e  -  ^{[G]T[H]->[G]}e{q  +  Aq}  }  (5.14) 

l e=l  e=l  J 

which,  in  a  condensed  form  is 


[kW  =  (RLt)  -  (RLF  (s-is) 

Many  numerical  examples  in  Chapter  3  and  4  prove  that  only  a  Newton-Raphson  iterative  solver  cannot 
obtain  the  entire  failure  solution  for  the  problems  with  damage,  especially  when  a  snap-back  appears  in  the 
load-deformation  curve. 

According  to  the  arc-length  method  proposed  in  (23,  24,  92],  an  unknown  loading  parameter  (A  +  d\)  is 
introduced  to  govern  the  load  increments.  Equation  (5.15)  is  modified  with  this  loading  parameter  as 

[K  W  =  (X3  +  <W){RLt}  -  {RfntF  (5.16) 

where  both  dX3  and  dq3  are  unknowns,  and  d\3  can  be  either  positive  or  negative.  The  orthogonality 
condition  (3.27)  is  chosen  to  be  the  constraint  equation  required  by  the  additional  unknown  d\3 . 
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5.2.3  Stability  conditions 


Following  the  stability  conditions  derived  for  displacement-based  and  stress-based  finite  element  approxima¬ 
tions  in  [5,  15,  111],  the  stability  conditions  of  the  stress-displacement  field  variational  problem  in  X-VCFEM 
are  stated  in  section  4.2.3.  They  are  positive  definite  [Hm]  and  [Hc],  unique  stress  interpolation  functions, 
and  non-zero  stress  parameters  for  all  non-rigid  body  displacement  fields.  The  two  conditions  can  be  sat¬ 
isfied  by  implementing  numerical  methods  in  section  4.5.  And  the  third  one  is  accomplished  by  choosing 
npm  >  riq  +  n™  +  n"  *  2  —  3  and  ngc  >  ncq-  3. 

5.3  Aspects  of  Numerical  Implementation 

5.3.1  Adaptive  criteria  for  cohesive  crack  growth 

A.  The  criterion  for  the  incremental  cohesive  crack  advance  into  matrix: 

The  static  deflection/ penetration  behavior  at  an  interface  has  been  the  subject  of  numerical  research  efforts 
in  the  past  years  and  many  significant  results  for  various  kinds  of  materials  have  been  obtained([3,  45, 
46,  63]).  The  fracture  toughness  ratio  of  the  interface  and  the  matrix  material  has  been  identified  as  the 
most  important  parameter  governing  the  crack  deflection/penetration  phenomenon.  Predicting  crack  growth 
requires  to  calculate  the  energy  release  rate,  G,  and  a  knowledge  of  the  surface  fracture  energy,  Gc.  In  this 
chapter,  we  denote  by  G'  and  G'c  the  energy  release  rate  and  the  critical  energy  release  rate  for  the  case  of 
growth  along  interfaces,  and  by  Gm  and  Grcn  the  corresponding  quantities  for  penetration  into  matrix. 

As  seen  in  previous  results,  stress  concentration  always  appears  in  the  matrices  around  fibers,  which  results 
in  cohesive  interfaces  between  fiber  and  matrix  becoming  weak  and  even  debonded.  Simultaneously,  damage 
at  the  interface  results  in  larger  concentrated  stress  fields  in  matrix.  Once  the  stress  in  matrix  reaches  some 
critical  value,  the  material  at  this  matrix  point  might  become  softening  and  damage  propagates  into  matrix 
from  the  interface.  All  points  with  critical  stresses  are  regarded  as  the  candidate  damage  position,  where  the 
criterion  is  necessary  for  selecting  the  crack  growth  direction,  along  the  interface  or  branching  into  matrix. 
The  candidate  positions  are  usually  chosen  from  the  Gaussian  integration  points  on  the  interface.  In  the 
program,  18  Gaussian  integration  points  are  distributed  between  any  two  consecutive  nodes  at  the  interface. 
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At  the  candidate  points,  the  criterion  for  assessing  the  crack  penetrating  into  matrix  is  defined  as 


Gm/Gi  >  G^/Gi 


(5.17) 


In  this  thesis,  the  energy  release  rate  is  calculated  based  on  cohesive  zone  models,  which  are  shown  in  figure 
5.2.  The  bilinear  model  in  figure  5.2  (a)  is  for  describing  the  damage  at  interface,  and  the  linear  model  in 
figure  5.2  (b)  is  for  the  matrix  cracking.  The  areas  of  the  shadow  regions  express  the  current  energy  release 
rates  Gl  and  Gm.  According  to  the  relation  between  the  cohesive  energy  (f>  for  complete  decohesion  and  the 
critical  energy  release  rate  Gc  in  equation  (4.56),  the  critical  release  rates  for  interface  and  matrix  are 
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(5.18) 


In  order  to  obtain  the  energy  release  rate  Gm,  the  effective  cohesive  traction  t  is  calculated  according  to 
stresses  ( axx ,  ayy  and  axy)  at  every  candidate  point.  Recalling  equations  (4.59-4.63)  in  Chapter  4,  effective 
cohesive  traction  t(ac),  the  cohesive  energy  <t>(ac)  and  the  energy  release  rate  Gm  are  obtained 


t{a)  = 


(<7xxsin2a  -  crxysin2a  +  cryycos2a)2  +  /?-2(—  -crxxsin2a  +  axycos2a  +  -cryysin2a)2 


(5.19) 


sm 

Gm  =  *(«c)  =  ~r~«aX2  ~  t(acf)  (5.20) 

max 

where  ac  is  the  angle  maximizing  the  cohesive  energy. 

The  current  energy  release  rates  for  interface,  G\  is  obtained 
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According  to  equations  (5.18,  5.20,  5.21)  and  inequality  (5.17),  the  damage  propagation  directions  are 
determined  at  the  candidate  points. 

B.  Direction  and  length  of  the  incremental  cohesive  crack  advance: 

Recalling  results  in  chapter  4,  the  direction  of  matrix  cracking  ac  is  obtained  at  the  damage  onset  points  as 


ac 


arctan 


(-<rxz+a-vv±y/(cTzx-al/v)2+  4<t^\ 
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(5.22) 


The  sign  in  equation  (5.22)  is  chosen  as  the  one  that  maximizes  the  cohesive  energy  <f>c  by  satisfying  the 
condition  in  equation  (4.62)  c. 

Upon  establishing  the  direction  of  incremental  cohesive  crack  growth  ac,  the  length  of  cohesive  zone  advance 
(A/)  should  be  estimated  in  the  crack  evolution  scheme  according  to  the  same  algorithm  shown  in  chapter 
4  as: 


Al=j-^±r\AB  I.  M 

<PA  ~  <PB 

where  B  is  a  point  close  to  A  in  the  direction  of  crack  propagation. 

5.3.2  Generation  of  [Gc] 

Once  damage  is  driven  from  interface  into  matrix,  two  node-pairs  (mi,  ni)  and  (m2,  n.2),  shown  in  figure 
5.3,  are  added  at  the  interface,  where  nodes  mj  and  m2  are  at  the  matrix  side  and  nodes  rii  and  ri2  are 
at  the  inclusion  side.  The  separation  between  mi  and  m2  describes  the  displacement  discontinuity  at  the 
crack  surface.  Since  crack  doesn’t  propagates  into  the  inclusion,  the  node-pair  (rii,  7x2)  merges  by  sharing 
the  same  displacement.  This  can  be  implemented  at  assembling  matrix  [Gc].  In  matrix  [Gc],  the  elements  in 
column  DOFn2  are  added  to  the  corresponding  elements  in  column  DOFn  1,  and  the  entire  column  DOFn 2 
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is  assigned  zero.  The  process  is  shown  in  equation  (5.24)  as 
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5.4  Numerical  Example 


(5.24) 


An  example  with  a  square  microstructure  containing  a  single  circular  fiber  with  a  debonding  interface  is 
considered  to  check  the  effectiveness  of  X-VCFEM  and  study  the  interaction  between  the  interface  and 
matrix  cohesive  cracking.  The  geometrical  dimensions  for  the  specimen  in  figure  5.4(a)  are  a=  20  mm, 
r=  5  mm.  The  material  parameters  for  matrix  and  fiber  are:  Young’s  modulus  Em  —  72  GPa,  Ej  = 
450  GPa,  and  Poisson  ratio  um  =  0.32,  vj  =  0.17,  where  subscript  (  )m  and  (•)/  denote  matrix  and 
fiber  respectively.  The  interface  uses  the  bilinear  cohesive  zone  model  with  the  properties  crlmax  =  0.04  GPa, 
Sc  —  0.001mm,  b\  =  0.02mm,  (P  =  0.707.  The  linear  cohesive  zone  model  is  used  to  describe  matrix  cracking 
with  parameters:  cr™ax  =  0.05  GPa,  5™  =  0.002mm,  and  f3m  =  1.  The  cohesive  parameters  are  chosen 
to  make  cr™ax  >  crlmax,  so  that  the  damage  starts  from  interface  instead  of  matrix.  Under  plane  strain 
conditions,  the  displacement  boundary  conditions  are  shown  in  figure  5.4(a).  The  whole  microstructure  is 
modeled  with  one  X-VCFEM  element,  consisting  of  16  nodes  on  the  cell  boundary  and  20  node  pairs  on  the 
interface  for  displacement  interpolation.  Before  damage  propagates  into  matrix,  the  stress  functions  in  this 
example  consist  of  102  terms  of  polynomial  functions  and  45  terms  of  reciprocal  functions.  After  the  cracks 
advance  into  matrix,  one  branch  function  and  16  wavelet  functions  are  added  into  the  stress  interpolation  for 
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each  crack.  Figures  5.4(b)  shows  the  contour  plots  of  the  microstructural  stress  ayy  together  with  evolved 
position  of  the  cracks  at  the  final  stage  of  loading.  The  growth  pattern  of  each  crack  can  be  observed  by 
comparing  with  its  initial  configuration  in  figures  5.4(a),  where  there  is  no  matrix  cracks. 

The  relation  of  the  propagation  of  multiple  cracks  to  the  interface  debonding  is  in  general  complicated. 
However,  several  observations  can  be  made  based  on  the  results  of  the  simulation  by  this  model. 

•  As  shown  in  figures  5.4(b),  the  lower  stress  at  point  A  implies  that  interface  there  becomes  weak  even 
debonded,  which  results  in  the  load  can  not  be  transfered  into  the  fiber  at  this  position  effectively.  The 
positions  with  concentrated  stress  bifurcate  from  point  A  and  move  to  left  and  right  sides  respectively 
along  the  interface,  which  might  drive  the  damage  into  matrix  from  the  interface.  The  same  thing 
happens  at  the  bottom  point  of  the  circular  fiber. 

•  Due  to  symmetry,  four  cracks  propagate  into  matrix  from  the  interface.  Largest  stresses  appear  at  tips 
of  cracks  and  the  stress  in  fiber  is  released.  In  this  example,  since  cracks  result  in  larger  concentrated 
stress  than  interfaces,  cracks  propagation  in  matrix  becomes  the  key  damage  phenomenon  in  following 
failure  process. 

•  The  evolved  crack  path  tends  to  align  in  a  direction  perpendicular  to  the  applied  load  direction,  which 
agrees  with  the  observation  in  chapter  4. 

5.5  Concluding  Remarks 

The  extended  Voronoi  cell  finite  element  model  is  improved  in  this  chapter  to  predict  the  damage  advanc¬ 
ing  into  matrix  and  study  the  interaction  between  interfacial  debonding  and  matrix  cracking.  Polynomial 
functions,  reciprocal  functions,  branch  functions  and  wavelet  functions  are  made  to  the  element  stress  in¬ 
terpolations  to  accurately  depict  the  stress  discontinuities  and  concentrations  at  interfaces  and  cracks.  The 
damage  in  interface  and  matrix  are  modeled  by  cohesive  models.  A  criterion  for  assessing  the  crack  pen¬ 
etrating  into  matrix  is  proposed,  which  is  based  on  the  energy  release  rate  and  cohesive  energy.  A  square 
specimen  containing  a  single  circular  fiber  with  a  debonding  interface  is  considered  to  check  the  effectiveness 
of  X-VCFEM  and  the  criterion. 
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The  damage  analysis  in  fiber-reinforced  composites  is  in  general  complicated.  X-VCFEM  is  easy  to  be  ex¬ 
tended  to  study  effects  of  material  properties  and  geometric  characterization,  such  as  clustering,  alignment, 
fiber  shape,  relative  sizes  etc.,  which  are  critical  to  the  failure  process  in  the  microstructure.  This  will  be 
explored  in  the  future  work. 
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Figure  5.2:  Cohesive  zone  models  for  calculating  energy  release  rates:  (a)  the  bilinear  law  for  interface 
debonding  and  (b)  the  linear  law  for  matrix  cracking. 


Chapter  6 

Concurrent  Multi-level  Model  for 
Damage  Evolution  in 
Microstructurally  Debonding 
Composites 

6.1  Introduction 

Analysis  of  composite  materials  with  microstructural  heterogeneities  is  conventionally  done  with  macroscopic 
properties  obtained  by  homogenizing  response  functions  in  the  representative  volume  element  (RVE)  from 
microscopic  analyses  at  smaller  length  scales.  While  these  “bottom-up”  homogenization  models  are  efficient 
and  can  reasonably  predict  macroscopic  or  averaged  behavior,  such  as  stiffness  or  strength,  they  have  limited 
predictive  capabilities  with  problems  involving  localization,  failure  or  instability.  Assumptions  of  macroscopic 
uniformity  and  RVE  periodicity,  the  two  basic  requirements  of  homogenization,  break  down  under  these 
circumstances.  The  uniformity  assumption  ceases  to  hold  in  critical  regions  of  high  local  solution  gradients, 
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such  as  near  free  edges,  interfaces,  material  discontinuities  or  evolving  damage.  RVE  periodicity,  on  the  other 
hand,  is  unrealistic  for  non-uniform  microstructures,  e.g.  in  the  presence  of  clustering  of  heterogeneities  or 
microscopic  damage.  Even  with  a  uniform  phase  distribution  in  the  microstructure,  the  evolution  of  localized 
stresses,  strains  or  damage  path  can  violate  the  periodicity  conditions.  Problems  like  this  have  been  effectively 
tackled  by  multi-scale  modeling  methods  e.g.  in  [81,  33,  50,  75,  74,  89,  83,  82,  101,  120,  99].  Multi-scale 
analyses  methods  can  be  broadly  classified  into  two  classes.  The  first  is  known  as  ’’hierarchical  models” 
[33,  50,  101,  99]  in  which  information  is  passed  from  lower  to  higher  scales,  usually  in  the  form  of  material 
properties.  The  hierarchical  homogenization  models  assume  periodic  representative  volume  elements  (RVE) 
in  the  microstructure  and  uniformity  of  macroscopic  field  variables.  The  second  class,  known  as  “concurrent 
methods”  [75,  74,  90,  83,  82,  120],  implement  sub-structuring  and  simultaneously  solve  different  models  at 
regions  with  different  resolutions  or  scales. 

The  two-way  coupling  of  scales  enabled  in  the  concurrent  methods  is  suitable  for  problems  involving 
localization,  damage  and  failure.  Macroscopic  analysis,  using  bottom-up  homogenization  in  regions  of  rel¬ 
atively  benign  deformation,  enhances  the  efficiency  of  the  computational  analysis.  As  a  matter  of  fact,  it 
would  be  impossible  to  analyze  large  structural  regions  without  the  advantage  of  a  continuum  model  based 
macroscopic  analysis.  On  the  other  hand,  the  top-down  localization  process  cascading  down  to  the  mi¬ 
crostructure  in  critical  regions  of  localized  damage  or  instability  for  pure  microscopic  analysis,  is  necessary 
for  accurately  predicting  the  damage  path.  These  microscopic  computations,  depicting  the  real  microstruc¬ 
ture  are  often  complex  and  computationally  prohibitive.  Hence,  a  concurrent  setting  makes  such  analyses 
feasible,  provided  the  ”zoom-in”  regions  are  kept  to  a  minimum.  The  adaptive  multi-level  models,  promoted 
in  [75,  74,  90,  83,  82,  120],  are  attempts  to  achieve  this  objective,  with  the  adaptivity  motivated  from  phys¬ 
ical  and  mathematical  perspectives.  However,  there  is  a  paucity  of  such  studies  in  the  literature  involving 
material  nonlinearity  and  evolving  microstructural  damage.  In  their  previous  studies,  Ghosh  and  coworkers 
have  proposed  adaptive  multi-level  analysis  using  the  microstructural  Voronoi  cell  FEM  model  for  modeling 
elastic-plastic  composites  with  particle  cracking  and  porosities  in  [89],  and  for  elastic  composites  with  free 
edges  and  stress  singularities  in  [83,  82], 

In  a  preceding  paper  [84],  the  authors  have  derived  and  computationally  modeled  an  anisotropic  con- 
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tinuum  damage  mechanics  (CDM)  model  for  unidirectional  fiber-reinforced  composites  undergoing  interfa¬ 
cial  debonding  from  by  using  homogenization  theory.  The  CDM  model  homogenizes  the  damage  incurred 
through  initiation  and  growth  of  interfacial  debonding  in  a  microstructural  RVE  with  nonuniform  distribu¬ 
tion  of  fibers.  Additionally,  arbitrary  loading  conditions  are  also  effectively  handled  by  this  model.  The 
present  paper  uses  this  CDM  model  of  [84]  in  an  adaptive  concurrent  multi-level  computational  model  to 
analyze  multi-scale  evolution  of  damage.  Damage  by  fiber-matrix  interface  debonding,  is  explicitly  modeled 
over  extended  microstructural  regions  at  critical  locations  [37,  58],  The  adaptive  model  addresses  issues  of 
efficiency  and  accuracy  through  considerations  of  physically-based  modeling  errors. 

The  adaptive  multi-level  model  consists  of  three  levels  of  hierarchy  viz.  level-0,  level- 1  and  level-2),  which 
evolve  in  sequence.  The  continuum  damage  model  developed  in  [84]  is  used  for  level-0  computations.  The 
level- 1  domain  is  used  as  a  ‘swing  region’  to  establish  criteria  for  switching  from  macroscopic  to  microscopic 
calculations.  Physical  criteria  involving  variables  at  the  macroscopic  and  microstructural  RVE  levels,  trigger 
switching  from  pure  macroscopic  to  pure  microscopic  calculations,  i.e.  the  level  -  0  — >  level  —  1  — » level  —  2. 
A  transition  layer  is  placed  between  the  level  -  1  and  microscopic  level  -  2  domains  for  smooth  transition 
from  one  scale  to  the  next.  All  computations  in  the  composite  microstructure  with  explicit  representations 
of  the  fiber  and  matrix  phases  are  done  with  the  Voronoi  cell  finite  element  model  or  VCFEM  [37,  58]. 
In  VCFEM,  debonding  at  the  fiber-matrix  interface  is  achieved  by  a  layer  of  cohesive  springs  [76].  Two 
numerical  examples  are  solved  in  this  paper  to  examine  the  effectiveness  of  the  multi-level  computational 
model  in  multi-scale  damage  analysis.  The  first  example  considers  a  small  region  of  a  fiber  matrix  composite 
microstructure  for  comparison  with  an  explicit  micromechanics  model.  The  second  set  of  problems  models 
a  double  lap  bonded  composite  joint  for  demonstrating  its  capability  in  handling  large  structural  problems. 

6.2  Levels  in  the  Multi-scale  Computational  Model 

The  multi-phase  composite  computational  domain  Clhet  is  adaptively  decomposed  into  a  set  of  non-intersecting 
open  subdomains,  belonging  to  levels-0,  -1  and  -2  with  different  algorithmic  treatments,  i.e.  flhet  = 
njoUniiUftalintr.  The  different  levels  of  computational  hierarchy,  in  the  order  of  sequence  of  emer- 
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gence,  are  depicted  in  figure  (6.1)  and  discussed  briefly  here. 


LEVEL  2 


(a)  (b) 

Figure  6.1:  Schematic  of  the  two-way  coupled  concurrent  multi-level  model:  (a)  a  representative  volume 
element  (RVE)  for  a  non-uniformly  distributed  composite  microstructure  generated  by  tessellating  the  lo¬ 
cal  microstructure,  (b)  the  top-down  multi-level  model  showing  components  of  concurrent  coupling,  viz. 
continuum  level-0 ,  level- 1  of  asymptotic  homogenization  and  lev  el- 2  of  micromechanical  analysis. 


6.2.1  Computational  Subdomain  Level-0  (O/o  ) 

This  level  corresponds  to  regions  where  continuum  constitutive  laws  can  be  used  in  macroscopic  analysis. 
Macroscopic  field  variables  like  stresses  and  strains  in  fi/o  are  relatively  uniform  and  there  is  no  strong  non- 
periodicity  in  the  microstructure.  Hence,  microscopic  ‘statistical’  periodicity  in  the  RVE  is  assumed  to  be 
valid  in  this  level.  Scale  effects  are  negligible  and  it  is  possible  to  derive  effective  constitutive  relations  by 
volume  averaging  the  RVE  response  with  imposed  periodicity  conditions,  in  the  limit  that  the  RVE  tends  to 
zero  volume.  This  is  generally  the  starting  level  in  the  multi-scale  analysis  model,  as  long  as  RVE’s  can  be 
identified  for  the  computational  domain.  Macroscopic  analysis  with  the  continuum  constitutive  models  in 
level-0,  reduce  the  computing  effort  by  several  orders  of  magnitude  in  comparison  with  models  that  require 
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complete  microscopic  analysis  . 

For  undamaged  microstructures  with  linear  elastic  or  elastic-plastic  phases,  homogenized  anisotropic 
constitutive  laws  have  been  developed  by  the  authors  in  [90,  36].  In  the  case  of  microstructures  with  randomly 
evolving  microcracks  causing  diffused  damage,  the  homogenized  material  behavior  is  best  represented  by  a 
continuum  damage  mechanics  (CDM)  law.  An  anisotropic  CDM  model  with  a  fourth  order  damage  tensor 
has  been  developed  from  rigorous  micro  mechanical  analyses  in  [84],  The  general  form  of  CDM  models 
[54]  introduce  a  fictitious  effective  stress  Ey  acting  on  an  effective  resisting  area  (A),  which  is  caused  by 
reduction  of  the  original  resisting  area  A  due  to  material  degradation  from  the  presence  of  microcracks  and 
stress  concentration  in  the  vicinity  of  cracks.  In  [84],  the  effective  stress  E,j  is  related  to  the  actual  Cauchy 
stress  Ey  through  the  fourth  order  damage  effect  tensor  Mx]ki  as 

Ey  =  A/ijfci(D)Efci  (6.1) 

where  Mijki  is  a  function  of  the  fourth  order  damage  tensor  D(=  D^kie-,  <8>  ej  ®  ®  ei).  The  hypothesis  of 

equivalent  elastic  energy  is  used  to  evaluate  Mijki  and  hence  establish  a  relation  between  the  damaged  and 
undamaged  stiffnesses  [22,  19,  118].  Equivalence  is  established  by  equating  the  elastic  energy  in  the  damaged 
state  to  that  in  a  hypothetical  undamaged  state  as 

W(X,D)  =  ±Xij(Eijklm-1Zk,  =  W(£,0)  =  (6.2) 

where  E  =  E^e;  ®  ej,  E°-kl  is  the  elastic  stiffness  tensor  in  the  undamaged  state  and  Eijki( D)  is  the  stiffness 
in  a  damaged  state.  From  equations  6.1  and  6.2,  the  relation  between  the  damaged  and  undamaged  stiffnesses 
is  established  as 

Eijki  =  (Mpqij )  ~ 1  E°qrs  {Mrski)- 1  (6.3) 

With  an  appropriate  assumption  of  a  function  for  M^ki,  equation  (6.3)  can  be  used  to  formulate  a  damage 
evolution  model  using  micromechanics  and  homogenization.  In  [84],  a  damage  evolution  surface  is  introduced 
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to  delineate  the  interface  between  damaged  and  undamaged  domains  in  the  strain  e-space  as 


F  =  i CijPijkieki  -  K{aWd)  =  0  (6.4) 

Here  Wd  corresponds  to  the  dissipation  of  the  strain  energy  density  due  to  stiffness  degradation  for  constant 
strain  without  an  external  work  supply.  Also  called  the  degrading  dissipation  energy  (see  [49]),  it  is  an 
internal  variable  denoting  the  current  state  of  damage,  and  is  expressed  as: 

Wd  —  j'  ~€ij£kidEijki  (6 ■  n) 

Pijki  is  a  symmetric  negative-definite  fourth  order  tensor  that  will  be  expressed  as  a  function  of  the  strain 
tensor  dj,  a  is  a  scaling  parameter  and  k  is  a  function  of  Wd.  Assuming  associativity  rule  in  the  stiffness 
space,  the  evolution  of  the  fourth  order  secant  stiffness  is  obtained  as 

dF 

Fijki  =  A— — j  r  =  A Pijki  (6-6) 

0\2eijekl) 

Pijki  (e)  corresponds  to  the  direction  of  the  rate  of  stiffness  degradation  tensor  Eijki  •  For  a  composite  material 
with  interfacial  debonding,  the  direction  of  rate  of  stiffness  degradation  varies  with  increasing  damage  and 
hence  Pijki  (e)  does  not  remain  a  constant  throughout  the  loading  process.  The  model  requires  the  evaluation 
of  k,  a  and  Pijki  in  equation  (6.4).  These  are  determined  from  the  results  of  micromechanical  simulations  of 
a  RVE  with  periodic  boundary  conditions.  The  function  n(Wd)  is  evaluated  for  a  reference  loading  path  and 
all  other  strain  paths  are  scaled  with  respect  to  this  reference.  Upon  determination  of  the  maximum  value 
Wd  for  a  reference  loading  condition,  the  value  of  a  for  any  strain  path  can  be  obtained  by  simple  scaling. 
To  account  for  the  variation  of  Pijki{e),  any  macroscopic  strain  evolution  path  is  discretized  into  a  finite  set 
of  points.  The  values  of  Pijki  are  explicitly  evaluated  at  these  points  from  RVE  based  simulations.  Values 
of  Pijki  for  any  arbitrary  macroscopic  strain  value  can  then  be  determined  by  interpolating  between  nodal 
values  using  shape  functions  of  a  3D  linear  hexahedral  element.  The  details  of  the  parameter  evaluation 
process  in  the  macroscopic  CDM  model  are  discussed  in  [84], 
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6.2.2  Computational  Subdomain  Level-1  (Qa) 

Level-1  is  an  intermediate  computational  subdomain,  introduced  as  a  swing  region  for  establishing  criteria 
for  switching  from  macroscopic  level-0  regions  to  level-2  regions  of  pure  microscopic  computations.  The 
switching  criteria  are  based  on  analyses  of  the  macroscopic  problem,  as  well  as  of  the  microstructural  RVE 
problem.  The  asymptotic  homogenization  theory  is  used  for  this  level  to  decouple  the  set  of  governing 
equations  into  a  set  of  (i)  homogenized  equations  representing  the  macroscopic  problem  corresponding  to  a 
length  scale  x,  and  (ii)  microscopic  equations  for  the  RVE  Y (x),  represented  by  a  length  scale  y.  Details  of 
the  decoupled  macro-  and  micro-equations  are  given  in  the  appendix  section  6.7.1. 

Gradients  of  important  field  variables  are  evaluated  from  macroscopic  analysis  to  assess  the  deviation 
of  macroscopic  uniformity.  Such  gradients  may  be  the  effect  of  strong  microscopic  non-homogeneity  in  the 
form  of  highly  localized  stresses  and  strains  or  damage.  The  RVE-based  microscopic  analysis,  on  the  other 
hand,  provides  effective  criteria  to  estimate  departure  from  periodicity  conditions,  especially  in  the  event  of 
evolving  microstructural  damage.  The  adaptation  criteria  for  level  transitions  are  discussed  in  section  6  4. 
Two  sets  of  finite  element  problems  are  solved  for  the  level- 1  subdomain  in  sequence,  viz., 

1.  Macroscopic  analysis:  Incremental  macroscopic  analysis  of  the  computational  domain  is  performed  using 
the  CDM  model  to  evaluate  macroscopic  variables  e.g.  stresses  and  strains  due  to  the  increments  in 
applied  loads. 

2.  Microstructural  RVE  analysis:  This  is  a  post-processing  operation  in  which  microstructural  analysis 
of  the  RVE  is  conducted  for  each  integration  point  of  the  macroscopic  elements.  The  strain  field  aj, 
obtained  from  macroscopic  analysis  with  the  CDM  model,  is  imposed  on  the  RVE  as  an  external  driver, 
together  with  periodic  boundary  conditions  on  the  boundary  of  the  RVE  as  shown  in  figure  (6.1)a. 
Microscopic  stresses  £7tJ,  strains  and  other  variables  are  computed  in  this  post-processing  stage  for 
each  RVE. 

Remark  1:  The  macroscopic  computations  of  level-0  and  level- 1  elements  are  performed  with  the  conventional 
displacement-based  finite  element  method,  while  all  microscopic  calculations  in  the  RVE  of  level-1  elements 
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are  performed  using  the  Voronoi  cell  FEM  [88,  37,  58]. 

Remark  2:  Computational  models  in  the  macroscopic  level-0  and  level- 1  subdomains  are  refined  adaptively  by 
selective  h-  or  h-p  strategies.  ‘Error’  and  convergence  criteria  for  this  refinement  have  been  discussed  in  [83]. 
Local  enrichment  through  successive  mesh  refinement  or  enhancement,  serves  a  dual  purpose  in  the  multi¬ 
level  computational  strategy.  The  first  goal  is  to  identify  regions  of  high  discretization  ‘error’  and  improve 
convergence  through  mesh  enhancement.  The  second  is  to  identify  regions  of  high  modeling  error  and  zoom 
in  on  these  regions  to  create  higher  resolution.  These  regions  are  generally  characterized  by  large  gradients 
and  localization  of  macroscopic  variables.  Element  refinement  in  these  regions  is  helpful  for  reducing  the 
length-scale  difference  between  macroscopic  elements  in  the  homogenized  domain  and  microscopic  regions 
with  explicit  representation  of  heterogeneities. 

6.2.3  Computational  Subdomain  Level-2  (fi/2  ) 

The  level-2  subdomain  of  pure  microscopic  analysis  emerges  from  level-1  elements  in  regions  characterized 
by  (a)  departure  from  macroscopic  uniformity,  e.g.  regions  of  localization  or  fracture,  and  (b)  significant 
microstructural  non-uniformities  manifested  by  e.g.  growth  of  localized  damage.  Prior  to  transition  to  level- 
2  elements,  a  high  spatial  resolution  is  reached  in  the  macroscopic  mesh,  resulting  in  small  elements,  by 
h-  or  hp-  refinement.  The  successive  refinement  process  stops  when  a  certain  element  size  is  achieved  and 
subsequently  the  model  changes  from  macroscopic  to  pure  microscopic.  A  scale  ratio  SR  is  chosen  a-priori 
to  ascertain  this  element  size.  Depending  on  the  choice  of  SR  —  SlZsi°{  of  local  »  the  microscopic 

model  in  any  given  level-2  element  can  encompass  large  portions  of  the  microstructure  with  many  discrete 
heterogeneities.  The  level-2  elements  are  constructed  by  filling  with  the  exact  microstructure  at  that  location, 
as  outlined  in  the  following  steps  and  shown  in  figure  6.2. 

•  Use  appropriate  adaptation  criteria  to  determine  if  a  level- 1  element  needs  to  switch  to  level- 2  element. 

•  Identify  a  region  in  the  microstructure  f2m;cro  that  is  located  in  the  same  region  as  the  level-2  element. 
Sl-micro  should  extend  beyond  the  element  boundary  by  approximately  two  fiber  lengths. 

•  Tessellate  the  local  microstructure  to  generate  a  mesh  of  Voronoi  cell  elements  as  shown  in  figure  (6.3). 
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•  Carve  out  the  microstructural  region  of  the  level-2  element  from  the  local  microstructure  fImicro.  This 
procedure  will  result  in  dissecting  some  of  the  fibers  on  the  boundary.  When  this  happens,  additional 
nodes  are  generated  on  the  Voronoi  cell  boundary  at  locations  where  the  fiber  surface  and  Voronoi 
cell  edges  intersect  the  boundary  of  the  level-2  element.  The  dissected  conjugate  pieces  of  a  fiber 
belonging  to  two  contiguous  level-2  elements  are  joined  together  when  the  two  contiguous  elements 
share  a  common  edge. 


<j> 

() 

O 

O 

© 

□  Level-0/1  Element  CO  Transition  element 

■  Level-2  element  G— ©  Special  interface  layer 

•  Level-0/ 1  nodes 

O  Level-0/1  nodes  at  the  transition  interface 

■  VCFEM  nodes  on  level-2/transition  boundary 
X  VCFEM  internal  nodes 

Y  Transition  element  nodes  at  the  interface 


(a)  (b) 

Figure  6.2:  (a)  Process  of  carving  out  level-2  element  microstructure  (b)  Interface  constraints  between  level- 
0 /level- 1  and  tr  elements 

Requirement  of  high-resolution  micromechanical  models  in  these  elements  entails  prohibitively  large  com¬ 
putations  using  conventional  finite  element  methods.  The  microstructure-based  Voronoi  cell  FEM  [37,  58,  88] 
is  particularly  effective  for  modeling  level-2  elements  because  of  its  efficiency  in  modeling  large  heterogeneous 
regions  [37,  58,  88,  89,  83].  Each  Voronoi  cell  with  embedded  heterogeneities  (particle,  fiber,  void,  crack  etc.) 
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represents  the  region  of  contiguity  for  the  heterogeneity,  and  is  treated  as  an  element  in  VCFEM.  VCFEM 
elements  can  be  considerably  larger  than  conventional  FEM  elements  and  incorporate  a  special  hybrid  FEM 
formulation.  Incorporation  of  known  functional  forms  from  analytical  micromechanics  substantially  enhances 
its  convergence.  A  schematic  diagram  of  Voronoi  cell  elements  is  shown  in  figure  (6.3).  A  high  level  of  ac¬ 
curacy  has  been  achieved  with  VCFEM  for  modeling  problems  with  microstructural  damage  by  particle 
cracking  [38]  and  fiber-matrix  interfacial  debonding  [37,  58],  For  debonding  simulation,  imperfect  interfaces 
are  represented  by  the  cohesive  zone  model  [76].  Displacement  degrees  of  freedom  on  the  fiber-matrix  in¬ 
terface  are  constrained  by  the  cohesive  zone  models  as  discussed  in  section  6.3.  VCFEM  has  been  shown 
to  be  significantly  more  efficient  than  commercial  displacement  based  FE  packages  for  modeling  complex 
microstructures  with  evolving  damage. 


Figure  6.3:  A  typical  level-2  element  containing  an  aggregate  of  microstructural  Voronoi  cell  elements  with 
relevant  notations. 
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6.2.4  Scale  Transition  Subdomain  (Qtr) 


The  interface  between  the  level-0  or  level-1  elements  and  the  level-2  elements  with  explicit  representation 
of  the  heterogeneous  microscopic  domain,  needs  a  special  treatment  to  facilitate  smooth  transition  of  scales 
across  the  element  boundaries.  A  layer  of  transition  elements  ( Etr  £  Otr)  is  sandwiched  between  these 
elements,  where  (Otr)  is  the  transition  subdomain  as  shown  in  figure  (6.2)b.  The  Etr  elements  are  essentially 
level-2  elements  with  compatibility  and  traction  continuity  constraints  imposed  at  the  interface  with  level- 
0/level-l  elements.  It  is  assumed  that  layers  of  Etr  elements  are  located  beyond  the  critical  hot-spots, 
at  which  homogenization  fails.  Hence,  the  homogenized  laws  are  sufficient  at  their  interfaces  with  level- 
l/level-0  elements.  A  weak  form  of  the  interface  displacement  continuity  is  incorporated  through  the  use  of 
Lagrange  multipliers  on  this  interface  [83,  82].  This  results  in  a  weak  satisfaction  of  the  interface  displacement 
compatibility  and  avoids  the  occurrence  of  spurious  forces  that  arise  if  the  displacements  are  strongly  coupled. 


6.3  Coupling  Different  Levels  in  the  Concurrent  Multi-Scale  Al¬ 
gorithm 

The  concurrent  multi-scale  analysis  requires  that  all  levels  be  coupled  for  simultaneously  solving  for  variables 
in  the  different  computational  subdomains.  Consequently,  the  global  stiffness  matrix  and  load  vectors  are 
derived  for  the  entire  computational  domain  (f Ihet  —  {Oio  U  On  U  fl(2  U  Htr}).  The  corresponding  domain 
boundary  is  delineated  as  =  {T/o  U  Ta  U  where  Tjo  =  dflio  D  I^et;  Fn  —  dtln  nTjej;  = 
80,12  O  r^et-  Let  Tint  =  80n  fl  80tr  delineate  the  boundary  between  the  level-1  and  transition  elements, 
where  the  displacement  continuity  is  satisfied  using  Lagrange  multipliers.  The  incremental  form  of  the 
equation  of  principle  of  virtual  work  equation  for  f l^et  at  the  end  of  an  increment,  can  be  written  as  the  sum 
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of  contributions  from  each  individual  domain,  as 
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/  (E tj  +  AEy)— i-  dfl  -  (U+  A U)5uf  dT 
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^  (ti  +  Ati)5uf  dr 
r12 
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+<5  [  (A‘r  +  AAf  )(ut  +  Av,  -  trf  -  Au‘r)dr  =  0 

J  rint 


(6.7) 


The  prefix  A  symbolizes  increments  of  the  respective  variables  in  the  incremental  solution  process.  The 
superscripts  10,  1 1,  12,  tr  correspond  to  association  with  the  respective  level,  while  the  (/)  sign  refers  to 
variables  that  could  belong  to  either  level.  E tJ-  are  the  components  homogenized  macroscopic  stresses  ob¬ 
tained  from  the  CDM  constitutive  model  for  O;o  and  fin.  The  applied  tractions  f*  are  at  traction  boundaries 
of  the  respective  domains.  The  boxed  parts  in  equation  (6.7)  correspond  to  contributions  from  level-2  and 
transition  computational  subdomains  that  are  generated  from  VCFEM  solutions  of  the  microstructural  re¬ 
gions.  Displacement  components  u[°,  uj1,  u[r  and  m[2  are  on  the  boundaries  of  elements  coinciding  with  the 
boundaries  of  the  flio,  fin,  fltr  and  fli i  subdomains.  An  intermediate  segment  rjnt  is  added  at  the  interface 
between  the  level- 1  and  tr  elements,  as  shown  in  figure  6.2.  On  these  segments,  displacement  components  Vi 
are  interpolated  with  any  order  polynomial  functions,  independent  of  the  interpolations  on  dfl ,0//il  or  dfltT . 
Even  for  highly  nonhomogeneous  displacements,  high  order  interpolations  on  the  intermediate  segment  are 
able  to  smoothen  the  transition  between  levels.  This  has  been  demonstrated  for  problems  without  damage 
through  numerical  examples  in  [83].  The  last  two  terms  in  equation  (6.7)  use  Lagrange  multipliers  to  fa¬ 
cilitate  incorporation  of  a  weak  form  of  the  interfacial  displacement  continuity  on  r;n(.  Xl0l/n  and  Xtr  are 
vector  columns  of  Lagrange  multipliers  belonging  to  domains  flio/n  and  fltr  respectively  at  rjn(.  The  Euler’s 
equations,  obtained  from  setting  the  coefficients  of  Svt,  <5Xll°''11  and  <5A(r  to  zero  respectively  in  the  principle 
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of  virtual  work  (6.7),  are 


A'0/n  +  AA-0/il  =  K-  +  Acrij),0/,1nj  =  -(A‘r  +  AA-r)  =  ~(ol3  +  A crij)trnj 

(«i  +  Awt),0/a  =  (ui  +  Au;)er  =  (uj  +  Avf)  (6.8) 

where  rij  is  the  unit  normal  vector  and  \f^n  and  A-r  correspond  to  the  interfacial  traction  components  on 
dflio/n  and  dfltr  respectively.  The  displacements  v,  and  the  Lagrange  multipliers  X[0^11  and  Xl2^tr  on  the 
intermediate  boundary  segment  are  interpolated  from  nodal  values  using  suitably  assumed  shape  functions 
as: 

{v}  =  [Lint}{qint}  ,  {Ai0/a}  =  [La‘o/u]{Aj0/ii}  >  {Atr}  =  [L^<r]{Atr}  (6.9) 

The  displacements  u-°  and  u[l  in  each  level- 0  and  level- 1  elements  are  interpolated  by  the  standard  or 
hierarchical  Legendre  polynomials  based  shape  functions  as: 

{u}'°  =  [Nio]{q<o}  =  [N[0  Np0] 

{u}I1  =  [Nil){q„}  =  [Nf1Ng] 

As  shown  in  figure  6.2,  the  generalized  displacements  in  the  level-0  and  level-1  elements  are  subdivided  into 
two  classes:  (i)  those  at  nodal  points,  which  interface  with  transition  elements,  and  (ii)  those  at  all  other 
nodes.  Generally,  only  level-1  elements  will  interface  with  transition  elements  because  of  the  sequence  of 
introduction  of  computational  levels.  The  generalized  displacements  q[0/,n  corresponds  to  the  nodal  degrees 
of  freedom  in  level-O/level-1  elements  at  the  interface  with  transition  elements,  while  q^/;i  correspond  to 
the  remaining  degrees  of  freedom  in  these  elements.  The  solution  of  the  algebraic  form  of  equation  (6.7)  is 
obtained  using  the  Newton  Raphson  iterative  solver.  Setting  up  the  tangent  stiffness  matrix  requires  consis¬ 
tent  linearization  by  taking  directional  derivative  of  equation  (6.7)  along  incremental  displacement  vectors 
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Au  and  Av,  and  the  Lagrange  multipliers  AA.  For  the  i— th  iteration  in  the  solution  of  the  incremental 


variables,  assembled  matrix  equations  derived  from  equation  (6.7)  has  the  following  structure. 


■ 

l 

i 

'  ' 

rv/0/M 

T {1,0 

^lo/n 

0 

0 

0 

P  io/n 

0 

Aqfo/u 

AFfo/u 

K-O.l 

^10/11 

* 

o  O 
20 

0 

0 

0 

0 

0 

Aqg/u 

AFpo/u 

0 

0 

I# 

Tfl.O 

0 

0 

Per 

Aq[r 

AF[r 

0 

0 

K?/ 

K0’0 
12/ tr 

0 

0 

0 

< 

Aq^/tr 

>  =  < 

AF°/tr 

0 

0 

0 

0 

0 

Q/o/u 

Qtr 

A  Qint 

AFinf 

PT 

rio/n 

0 

0 

0 

Q  lo/n 

0 

0 

AAio/u 

AFA(0/n 

0 

0 

Pr 

*  tr 

0 

Q  Tr 

0 

0 

AAtr 

J 

AFA12/tr 
.  > 

As  explained  before,  the  superscript  I  represents  quantities  associated  with  nodal  points  at  the  interface  with 
transition  elements  while  superscript  O  indicate  association  with  nodes  at  other  regions.  The  two  notations 
in  the  superscript  separated  by  comma,  represents  the  node  coupling  effect.  For  example,  the  superscript 
I,  O  corresponds  to  the  coupling  between  the  non-interface  and  interface  nodes.  The  stiffness  submatrices 
[K/o/a]  and  sub-vector  (F;o/n }  correspond  to  those  for  the  level-0  and  level-1  elements  and  are  expressed  as 


{Kl0/ll)man(3  —  / 

J  n. 


dNa  <9£mn  ON0 
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(6.12) 


The  subscripts  (m,  n )  correspond  to  the  degrees  of  freedom  while  (a,  /?)  correspond  to  the  node  numbers  in 
the  element.  These  matrices  and  vectors  are  further  divided  based  on  the  classification  of  the  I  and  O  nodes. 
The  coupling  between  the  level- 0 /level- 1  and  tr  elements  is  achieved  through  the  [P]  and  [Q]  matrices,  which 
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may  be  expressed  as 


( PlO/l\)mcmp  —  /  ^mo(^'Xl0/,1)n/3^' 

Jrint 

(Ptr)manp  —  ~  N^Q(LAi2/tr-)n/3dr 

'Tint 

(QlO/ll^macnft  =  I  (IJjnj)TnQr(Ij^io/ii  )n/?dF 

Jrint 

(Q  tr)manf 3  —  /  (IJjrjt)mor(I'v)n/}dT 
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(6.13) 


Contributions  to  the  load  vector  {F}  due  to  coupling  between  level- 0/level- 1  and  tr  elements  are  given  as 


(AFinf)mQ  =  -  /  (LltUX10'11  +  A\l0'n)mdr  -  [  (Lfnt)a(\l2/tr  +  A\n'tr)mdT 

int  J  rjnt 

(^^AIO/il)mo  ==  /  )a{^m  “f"  At/m  (^{0/il)m  A(u/o/Zl  )m 

(Ai^\|2/tr)ma  ~  I  (■^/Atr  (^Z2/Zr)m  —  A(uj2/tr)m}^r' 

''Tint 


(6.14) 


Finally,  the  stiffness  [KI2/tr]  and  the  load  vector  {F|2/tr}  for  level-2  and  tr  elements  are  obtained  by  VCFEM 
calculations  followed  by  static  condensation  to  represent  the  virtual  work  in  terms  of  the  boundary  terms 
only. 


6.3.1  Modified  Voronoi  Cell  FEM  Formulation  for  a  RVE  in  Level- 1  Elements 

Details  of  the  Voronoi  Cell  FEM  are  provided  in  [88,  37,  58]  and  have  been  summarized  in  the  appendix  section 
6.7.2.  As  discussed  in  section  6.2.2,  the  post-processing  phase  for  level-1  elements  require  the  evaluation  of 
different  variables  in  the  RVE  from  known  values  of  macroscopic  strains.  A  small  variant  of  the  formulation  in 
equation  (6.40)  enables  this  execution.  The  energy  functional  for  a  RVE  ( F )  with  F-periodic  displacements 
and  F-anti-periodic  tractions  on  the  boundary,  and  imposed  macroscopic  strain  (e^  +  Ae,j),  may  be  written 
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as 
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(6.15) 


(6.16) 


The  boxed  term  corresponds  to  the  additional  energy  due  to  the  imposed  macroscopic  strain  field  on  the 
RVE  region  Y.  The  Euler- Lagrange  equations  corresponding  to  this  energy  functional  are: 


<Mx>y)  +  Acij(x,y)  =  Sijki{<Tij  +  A  Oij)  =  (eij(x)  +  Aey(x)) 
,  lfd(u.(y)  +  Auj(y))  ,  d(uj(y)  +  Auj(y)), 


Ui  is  T-periodic  and  471)  *s  ^-anti-periodic  on  dYe 


(6.17) 

(6.18) 


The  corresponding  weak  form  of  the  element  kinematic  relation  is  written  in  a 


matrix  equation  form  as 
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(6.19) 
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or  in  a  condensed  form 


[He]{/3  +  A/3}  =  [Ge]{q  +  Aq}  -  {R f  }  (6.20) 

This  relation  is  then  substituted  in  equation  (6.44)  for  obtaining  the  RVE  based  solutions.  It  should  be 
noted  that  displacement  periodicity  is  imposed  on  the  RVE  boundary  for  solving  this  problem. 

6.4  Criteria  for  Adaptive  Mesh  Refinement  and  Level  Transitions 

In  the  application  of  the  multi-level  model,  the  following  criteria  are  used  for  mesh-refinement  and  level 
transitions  due  to  discretization  and  modeling  error  respectively.  Many  of  these  adaptation  criteria  are 
based  on  the  physics  of  the  problem  in  consideration,  since  rigorous  mathematical  error  bounds  are  scarce 
(or  even  non-existent)  for  these  nonlinear  problems.  Consequently  they  are  nonunique  and  other  indicators 
may  be  used  if  appropriate. 

6.4.1  Refinement  of  Level-0  and  Level-1  Meshes  by  h-Adaptation 

The  computational  models  in  level-0  and  level-1  subdomains  are  enriched  by  h-adaptation  based  mesh  re¬ 
finement  to  reduce  discretization  ‘error’.  The  h— adaptation  procedure  subdivides  candidate  macroscopic 
elements  into  smaller  elements  to  reduce  a  suitably  chosen  error.  It  is  necessary  to  impose  boundary  dis¬ 
placement  compatibility  constraint  conditions  between  contiguous  divided  and  undivided  elements  in  this 
method  [90].  This  local  mesh  enrichment  is  intended  to  reduce  discretization  error  and  to  identify  regions 
of  modeling  error  by  zooming  in  on  localization  regions  with  evolving  gradients.  For  CDM  based  evolving 
problems,  an  adaptation  criterion  is  formulated  in  this  paper  in  terms  of  the  jump  in  traction  across  adjacent 
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element  boundaries  that  signifies  local  stress  gradients.  The  condition  is  stated  as: 


Refine  element  ‘e\  if  the  traction  jump  error  across  the  element 
satisfies  the  condition:  E\]  >  C\  *  E^vg, 


where 


ptj  rE£?(iff)2 
av 9  1 


)1/2  and  (E^)l  = 


fan.  m)?  +  [[Ty\}2)ddn 

fdO'dM 


(6.21) 


Here  NE  is  the  total  number  of  level-0  and  level- 1  elements  in  the  entire  computational  domain,  Tx,Ty  are 
the  components  of  element  boundary  tractions  in  the  x  and  y  directions  and  [[.]]  is  the  jump  operator  across 
element  boundary  dQ.e.  A  factor  C i  (<  1)  has  been  chosen  from  numerical  experiments. 


6.4.2  Criteria  for  Switching  from  Level-0  to  Level-1  Elements 

Level-0  to  level- 1  element  transition  takes  place  according  to  criteria  signaling  departure  from  conditions  of 
the  homogenizability  that  are  based  on  macroscopic  variables  in  the  continuum  model  of  level-0  elements. 
The  degrading  dissipation  energy  Wd  in  the  CDM  model  is  a  strong  indicator  of  localized  damage  evolution. 
Consequently,  a  criterion  is  formulated  as: 


Switch  element  ’e’  from  level-0  to  level-1  if  : 

Efe  *  (Wd)e  >  C2  *  EtL  *  (Wd)mar  (6.22) 

where  Efde  is  the  norm  of  the  local  gradient  of  the  degrading  dissipation  energy  (Wd)e,  expressed  as: 

Efe  =  »/(ffi^)  2  +  (^p^)2  (6.23) 

E^x  *s  the  maximum  value  of  E§de  for  all  elements  and  C2  (<  1)  is  a  prescribed  factor.  The  criterion  (6.22) 
is  helpful  for  seeking  out  regions  with  high  gradients  of  Wd  in  regions  of  high  Wd  itself.  In  a  previous  paper 
by  the  authors  [83],  the  gradient  E-!dc  was  expressed  in  terms  of  the  maximum  difference  in  the  damage  for 
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neighboring  elements  as  Efde  —  Max|(Wd)e  —  (Vi^)  adjacent  I-  A  more  accurate  definition  of  the  local  gradient 
is  adopted  in  the  present  work,  using  the  Zienkiewicz-Zhu  (ZZ)  gradient  patch  recovery  method  [119].  In 
this  method,  interpolation  of  Wd  is  assumed  in  the  form  of  a  polynomial  over  a  patch  of  elements  adjoining 
a  nodal  point  in  a  level-0  element.  The  least  square  minimization  process  leads  to  the  local  matrix  equation 
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(6.24) 


where  [Ne(xi,  2:2)]  is  a  matrix  containing  polynomial  interpolation  terms  and  ne  is  the  number  of  elements 
in  the  patch.  The  equation  (6.24)  is  solved  for  the  coefficients  {a}.  The  gradients  of  Wd  in  each  element  are 
calculated  from  the  nodal  values  using  element  shape  functions  as 
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6.4.3  Criteria  for  Switching  from  Level-1  to  Level-2  Elements 

For  elements  in  which  macroscopic  nonuniformity  has  been  established  according  to  equation  (6.22),  de¬ 
parture  from  RVE  periodicity  is  taken  as  an  indicator  for  activating  a  switch  from  level-1  to  level-2.  The 
switching  criterion  is  developed  in  terms  of  evolving  variables,  e.g.  the  averaged  strain  at  the  fiber-matrix 
interface  in  the  local  microstructural  RVE.  The  averaged  strain  is  stated  as: 
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where  the  integral  is  evaluated  over  all  the  fiber-matrix  interfaces  in  the  RVE.  The  jump  in  displacement 
across  the  fiber-matrix  interface  with  a  normal  rii  is  denoted  by  [uj].  For  perfect  interfaces  [uj]  will  be  zero. 
Thus,  Dij  corresponds  to  the  contributions  to  macroscopic  strain  due  to  damage  only,  and  Dij  =  0  in  the 
absence  of  damage.  Departure  from  periodicity  will  result  in  a  significantly  different  averaged  strain  Dij 
in  response  to  different  conditions  on  the  boundary  of  the  microstructural  region.  For  example,  let  D ®ji2 
correspond  to  the  solution  of  a  boundary  value  problem  of  the  local  microstructure  included  in  a  level-2 
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element  (see  figure  6.2)  subject  to  boundary  displacements  that  have  been  obtained  from  the  macroscopic 
level-0/1  analysis.  The  scale  of  the  microstructure  is  relevant  in  this  analysis  since  periodicity  is  not  imposed 
on  the  boundary.  On  the  other  hand,  let  D E  be  from  the  solution  of  a  boundary  value  problem  of  the 
local  RVE  with  imposed  macroscopic  strains  and  subjected  to  periodic  boundary  displacements  constraints. 
The  difference  in  these  two  strains  for  a  level- 1  element  e  may  be  quantified  as 
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For  evaluating  D^j2  in  a  given  step  of  the  incremental  solution,  only  the  increments  in  the  present  step 
are  calculated  by  the  level-1  macroscopic  displacement  boundary  conditions.  It  is  assumed  that  the  R.VE- 
based  solution  is  valid  all  the  way  upto  (but  excluding)  the  present  step.  The  departure  from  periodicity  is 
measured  in  terms  of  the  difference  in  averaged  strains  Edper,  and  hence  the  criterion, 


Switch  element  V  from  level-1  to  level-2  if : 

Edr r  >  C3DZxE  (6.28) 

where  DZx  maximum  value  of  \D^RVE\  for  all  the  level-1  elements  in  the  computational  domain. 

Remark:  Once  the  regions  of  level-2  and  transition  elements  have  been  identified,  it  is  important  to  update 
the  local  micromechanical  states  of  stress,  strain  and  damage  to  the  current  state.  This  step  should  precede 
the  coupled  concurrent  analysis.  For  this  analysis,  the  history  of  the  macroscopic  displacement  solution  on 
the  levcl-O/level-1  element  boundary  prior  to  the  switch  is  utilized.  The  local  micromechanical  (VCFEM) 
boundary  value  problem  for  the  level-2  element  is  incrementally  solved  from  the  beginning  to  obtain  the 
history  of  stresses,  strains  and  damage  in  the  microstructure  from  the  macroscopic  boundary  displacement 
history. 
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6.5  Numerical  Examples  with  the  Adaptive  Multilevel  Model 


Two  sets  of  numerical  examples  are  solved  to  study  the  effectiveness  of  the  multi-level  computational  model 
for  composite  materials. 

6.5.1  Multi-level  Model  vs.  Micromechanical  Analysis 

This  example  is  aimed  at  understanding  the  effectiveness  of  the  multi-level  model  in  analyzing  a  nonuniform 
composite  microstructure  by  comparing  its  predictions  with  those  by  pure  micromechanical  analysis.  It  is 
computationally  intensive  to  conduct  pure  micromechanical  analysis  with  evolving  damage  for  very  large  mi- 
crostructural  regions.  Consequently  a  computational  domain  with  a  small  population  of  fibers,  as  shown  in 
the  optical  micrograph  of  figure  6.4(a),  is  considered.  The  micrograph  is  for  a  polymer  matrix  composite  with 
a  random  dispersion  of  uniaxial  fibers.  The  dimensions  of  the  micrograph  analyzed  are  100/im  x  70.09 /im, 
containing  264  circular  fibers  of  diameter  1.645 fim  each,  corresponding  to  a  volume  fraction  of  32%.  Though 
the  domain  may  not  be  adequate  for  a  clear  separation  between  continuum  and  micromechanical  regions 
(since  relatively  large  regions  are  needed  to  materialize  the  RVE),  the  results  of  this  example  are  enough  to 
show  the  effectiveness  of  the  overall  framework. 

The  optical  micrograph  is  mapped  onto  a  simulated  microstructure  with  circular  fibers  that  is  tessellated 
into  a  mesh  of  264  Voronoi  cell  elements,  as  shown  in  figure  6.4(b).  The  constituent  materials  in  the  composite 
system  are  an  epoxy  resin  matrix,  stainless  steel  reinforcing  fibers  and  a  very  thin  film  of  freekote  (<  0.1  /im.) 
at  the  fiber-matrix  interface.  The  freekote  imparts  a  weak  strength  to  the  steel-epoxy  interface,  which  allows 
a  stable  growth  of  the  debond  crack  for  experimental  observation.  The  experimental  methods  of  material  and 
interface  characterization  have  been  discussed  in  [37].  Both  the  matrix  and  fiber  materials  are  characterized 
by  isotropic  elasticity  properties.  The  matrix  material  has  a  Young’s  modulus,  Eepoxy  =  4.6  GPa  and 
Poisson’s  ratio,  uepoxy  —  0.4,  while  the  fiber  material  has  a  Young’s  modulus,  Esteei  =  210  GPa  and 
Poisson’s  ratio,  i/s(ee;  =  0.3.  A  bilinear  cohesive  law  described  in  [58,  76]  is  used  in  this  analysis  for  modeling 
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(a)  (b) 

Figure  6.4:  (a)  Optical  micrograph  of  a  steel  fiber-epoxy  matrix  composite  with  264  fibers  (b)  the  simulated 
computational  model  with  a  Voronoi  cell  mesh 


the  fiber-matrix  interface.  In  this  model,  the  normal  and  tangential  tractions  are  given  as 
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where  t  is  a  bilinear  function  of  the  interfacial  separation  as 


(6.29) 


t  = 


Oc 


V<5  <  6C 


6-6,  _ 

(  &c-6,  °rnax 


V<5  >  5C 


(6.30) 


The  unloading  behavior  in  the  hardening  region  is  linear  following  the  loading  path.  In  the  softening  region, 
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the  unloading  proceeds  along  a  different  linear  path  from  the  current  position  to  the  origin  with  a  reduced 
stiffness,  for  which  the  t  —  S  relation  is 


t  = 


■'max  Umax 


Sc  Sjnax  ^  Se  and  5  S ^ 


(6.31) 


It  is  expected  that  the  degrading  dissipation  energy  Wj  in  the  macroscopic  CDM  model  depends  on  the 
cohesive  parameters  in  the  microstructural  debonding  model.  A  square  RVE  with  a  single  circular  fiber  is 
simulated  for  interfacial  debonding  with  three  different  sets  of  cohesive  parameters,  as  shown  in  the  inset 
of  figure  6.5.  The  cohesive  energies  are  the  same  for  all  cases.  However  in  one  case,  the  critical  separation 
length  Se  is  increased  while  in  the  other,  the  corresponding  peak  stress  omax  is  reduced.  The  figure  6.5  infers 
that  while  Se  has  a  small  influence  on  Wj,  the  effect  of  amax  is  certainly  significant,  at  least  in  the  early 
stages  of  straining. 


The  cohesive  parameters  used  in  this  paper  are:  amax  —  0.005  GPa,  6C  =  5.1  x  10~5  m  and  Se  — 
3.1  x  10~4  m.  The  microstructure  is  loaded  in  tension  in  the  horizontal  direction  with  a  displacement  of 
0.1  /xm  that  is  uniformly  increased  in  20  equal  increments,  corresponding  to  a  total  strain  of  eu  =  0.1%.  The 
displacement  boundary  condition  is  imposed  along  the  right  edge,  as  shown  in  figure  6.4(b). 

Micromechanical  analysis  by  VCFEM 

The  pure  micromechanical  VCFEM  solution  using  the  mesh  of  figure  6.4(b)  has  been  presented  in  [58]  and 
are  used  here  as  reference  solutions  for  the  multi-scale  simulation.  Figure  6.6(a)  shows  the  contour  plot  of 
microscopic  stress  crxx  at  the  final  step  of  the  micromechanical  simulation  with  a  depiction  of  interfacial 
debonding.  The  right  side  of  the  microstructure  shows  significant  concentrated  damage  with  this  load.  The 
debonding  initiates  at  the  top  and  percolates  to  the  bottom  of  the  microstructure  along  a  narrow  band. 

Multi-Scale  Analysis  with  the  Multi-Level  Model 

Multi-scale  analysis  is  performed  by  the  concurrent  multi-level  computational  model  and  the  results  are 
compared  with  those  from  the  micromechanical  VCFEM  analysis.  For  the  multi-level  model,  the  entire  corn- 
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Figure  6.5:  The  degrading  dissipation  energy  Wd  as  a  function  of  strain,  evaluated  for  different  cohesive  zone 
parameters  in  the  bilinear  cohesive  law 

putational  region  of  264  fibers  is  first  divided  into  9  macroscopic  finite  elements  as  shown  in  figure  6.7(a).  For 
evaluating  the  homogenized  constitutive  properties  for  each  of  element,  statistically  equivalent  representative 
volume  element  or  SERVE  for  the  microstructure  underlying  each  macroscopic  element  is  first  identified. 
Various  statistical  methods  have  been  used  to  determine  the  size  scale  of  the  RVE  and  the  number  of  inclu¬ 
sions  contained  in  it  [79,  39,  94,  1 06] .  Rigorous  methods  of  evaluating  statistically  equivalent  representative 
volume  elements  by  a  combination  of  statistical  methods  and  micromechanical  analyses  have  been  conducted 
by  the  first  author  in  [91,  97).  However,  since  the  number  of  fibers  in  the  micrograph  is  limited  in  this  ex¬ 
ercise,  a  simpler  assumption  is  made.  The  SERVE  for  each  element  is  assumed  to  consist  of  all  the  fibers 
belonging  to  that  element.  For  example,  to  generate  the  SERVE  for  an  element  window  in  the  micrograph 
of  figure  6.4(b),  all  fibers  whose  centers  are  located  within  this  window  are  first  identified  as  constituents 
of  the  RVE.  This  is  shown  by  the  aggregate  of  black  fibers  in  figure  6.8(a).  The  homogenization  method, 
discussed  in  sections  6.2.2  and  6.3,  requires  a  periodic  distribution  of  the  RVE  and  this  is  achieved  by  locally 
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repeating  the  arrangement  of  fibers  in  both  the  xi  and  xi  directions  for  a  period  length  in  figure  6.8(a).  This 
means  that  for  each  fiber  identified  in  the  element,  at  (x\,  X2),  four  identical  fibers  are  placed  at  the  locations 
(xi  ±  Xi,  X2),  (xi,X2  ±  X2)  where  (Xi,  X2)  are  periods  in  the  two  directions.  The  period  lengths  X\,  X2  are 
selected  such  that  the  volume  fraction  of  RVE  matches  that  of  the  local  microstructure.  Finally,  the  domain 
is  tessellated  into  a  network  of  Voronoi  cells  as  shown  in  figure  6.8(b)  Tessellation  provides  a  natural  way 
of  creating  periodic  SERVE  boundary.  For  non-uniform  fiber  arrangements,  the  SERVE  boundary  consists 
of  non-straight  line  edges.  The  nodes  on  this  SERVE  boundary  are  periodic,  i.e.  for  every  boundary  node 
a  periodic  pair  can  be  identified  on  the  boundary  at  a  distance  of  one  period  along  each  of  the  coordinate 
directions.  In  figure  6.8(b),  the  node  pairs  are  identified  as  A  A,  BB  etc.  The  number  of  fibers  and  their 
distribution  in  the  SERVE  of  each  macroscopic  element  is  shown  in  figure  6.7(a). 

Since  the  number  of  elements  in  this  exercise  is  very  small  (only  9),  level-0  simulations  with  the  CDM 
model  is  bypassed  in  the  multi-level  analysis.  All  elements  are  level- 1  at  the  start  of  the  multi-level  simulation. 
Switch  to  level-2  elements  is  made  in  accordance  with  equations  (6.27)  and  (6.28)  with  C2  =  0.2.  However 
the  D*j2  -  DtjRVE  terms  for  each  element  in  equation  (6.27)  are  replaced  by  the  difference  in  RVE  based 
averaged  strains  between  adjacent  elements  D'-j ' RV  E  -  D\2,RVE .  Also,  as  opposed  to  an  entire  macroscopic 
element,  a  single  layer  of  transition  Voronoi  cell  elements  is  included  between  the  level- 1  and  level-2  elements. 
In  figure  6.7(b)  the  Voronoi  elements  containing  the  grey  fibers  constitute  the  transition  layer,  while  those 
containing  the  black  fibers  belong  to  level-2.  An  interface  segment  T int  is  inserted  between  the  transition  and 
level-1  elements  at  a  distance  LtT/i2  from  the  right  edge.  Convergence  properties  of  the  multi-level  model 
are  studied  by  considering  two  cases  with  L‘)(12  =  0.35  and  Ltr/12  =  0.45.  This  is  achieved  by  changing  the 
size  of  the  initial  level-1  elements. 

As  depicted  in  figure  6.7(b),  only  three  elements  (3,  6  and  9)  at  the  right  side  of  the  initial  mesh  switch 
from  level-1  to  level-2.  A  comparison  of  results  by  (a)  VCFEM  based  micromechanical  analyses  (all  level-2 
elements)  ,  (b)  homogenization  based  macroscopic  analysis  (all  level-1  elements),  and  (c)  concurrent  multi¬ 
level  analysis  ( level-1  and  level-2  elements)  for  L‘r/‘2  =  0.35  and  0.45  is  made.  Contour  plots  of  crn  (GPa) 
showing  interfacial  debonding  at  the  end  of  the  simulation  are  shown  for  the  concurrent  multi-scale  analysis 
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in  figure  6.6(b,c).  The  discrepancy  in  the  damage  path  predicted  by  the  micromechanical  analysis  and  the 
multi-level  analysis  reduces  sharply  with  increasing  value.  This  can  be  attributed  to  the  fact,  that 

the  damage  path  is  very  sensitive  to  the  macro-micro  interface  conditions.  Since  the  sample  size  is  small 
and  there  is  no  real  periodicity  in  the  microstructure,  the  proximity  of  the  level-1  boundary  to  the  damage 
localization  zone  alters  the  local  boundary  conditions.  However  as  this  distance  is  increased,  the  microscopic 
stress  distribution,  debonding  pattern  and  damage  zone  replicates  the  real  event  observed  in  micromechanical 
analysis.  The  distribution  of  the  micromechanical  stresses  a n ,  generated  by  pure  micromechanical  and  multi¬ 
level  analyses,  are  plotted  along  a  line  through  the  middle  of  micrograph  in  figure  6.9.  The  micromechanical 
stresses  show  only  minor  oscillations  about  an  averaged  value  of  the  0.005  GPa  in  the  region  to  the  left 
of  the  level- 1 -level-2  interface.  In  the  region  to  the  right,  where  damage  is  predominant,  there  is  clearly  a 
convergence  of  the  stresses  with  increasing  Ltr/in  value. 


The  macroscopic  or  averaged  stress-strain  response  for  element  1  (always  level- 1)  and  element  9  (changes 
levels)  are  plotted  in  figures  6.10.  For  the  micromechanical  problems  with  debonding,  the  volume  averaged 
stresses  and  strains  are  evaluated  by  averaging  the  local  fields  over  the  microscopic  domain  as: 


Ejj  =  —  [  Oij(x\,xi)d£l  and  =  -=-  f  eij(xi,X2)d£l  -  Di 

“  Jn  “  Jn 


(6.32) 


where  Dij  is  the  strain  jump  defined  in  equation  (6.26).  The  results  for  all  the  models  are  in  good  agreement 
for  the  element  1,  where  there  is  no  significant  microstructural  damage.  The  small  difference  is  due  to  the 
periodicity  constraints  imposed  on  the  microstructure.  Also  there  is  a  difference  between  the  results  of  case 
1:  '2' 12  =  0.35  and  case  2:  =  0.45,  due  to  the  interface  conditions  at  rinj.  However,  as  is  expected 

the  results  are  quite  different  for  element  9,  where  significant  damage  is  observed  in  figure  6.6.  The  level-1 
analysis  shows  significant  deviation  from  the  micromechanical  analysis  due  to  imposed  periodicity  in  the 
damage  zone.  Once  again,  the  results  improve  significantly  with  increasing  Ltrln  ratio. 
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6.5.2  A  Composite  Double  Lap  Joint  with  Microstructural  Debonding 

Adhesive  bonded  joints  are  considered  as  preferred  alternatives  to  fasteners  for  joining  structural  components 
due  to  their  light  weight.  However,  adhesively  bonded  structures  consisting  of  different  materials,  can  induce 
high  stresses  near  the  interface  leading  to  failure  initiation  by  interfacial  debonding.  A  double-lap  bonded 
joint  with  boron-epoxy  composites  as  adherents,  is  analyzed  in  this  example.  An  adhesive  shown  as  ABCD 
in  figure  6.11(a)  is  used  to  bond  the  two  composite  materials.  Only  a  quarter  of  the  joint  is  modeled  from 
considerations  of  symmetry  in  boundary  and  loading  conditions.  For  boundary  conditions,  the  displacement 
component  u i  is  set  to  zero  along  the  face  x2  =  0  implying  symmetry  about  the  X\  axis.  The  displacement 
components  U\  and  u2  along  the  face  sq  =  8h  are  set  to  zero  corresponding  to  a  fixed  edge.  A  tensile 
displacement  rq  is  applied  on  the  face  of  the  lower  ply  at  aq  =  0.  Both  plies  above  and  below  the  adhesive 
are  made  of  unidirectional  boron  fiber-  epoxy  matrix  composite  materials.  The  fibers  are  uniformly  arranged 
in  a  square  array  in  the  microstructure,  implying  a  square  unit  cell  with  a  single  circular  fiber.  The  epoxy 
matrix  has  a  Youngs  modulus  E  =  4.6  GPa  and  Poisson’s  ratio  u  =  0.4,  while  boron  fibers  have  a  Youngs 
modulus  E  =  210  GPa  and  Poisson’s  ratio  u  =  0.3.  The  material  properties  of  the  isotropic  adhesive  are: 
Young’s  modulus  E  =  3.45  GPa  and  Poisson’s  ratio  v  =  0.34.  The  bilinear  cohesive  law  parameters  for  the 
matrix-fiber  interface  are:  amax  =  0.02  GPa,  5C  =  5.0  x  10-5  m  and  Se  —  20.0  x  10~4  m. 

Multi-level  analysis  for  model  with  450  fibers 

In  this  model,  the  top  ply  above  the  adhesive  consists  of  10  rows  of  fiber,  while  the  bottom  row  consists  of 
5  rows  resulting  in  a  total  of  450  fibers.  The  microstructural  volume  fraction  of  fibers  is  V/  =  20%.  The 
applied  displacement  on  the  face  at  Xi  =  0,  is  uniformly  increased  from  zero  to  u\  =  1.2  x  10~3/i  in  15 
uniform  increments.  The  number  of  fibers  is  kept  low  in  this  example,  such  that  a  micromechanical  analysis 
can  be  easily  done  for  this  example  with  a  mesh  of  450  Voronoi  elements,  each  of  which  is  a  square  unit  cell. 
The  micromechanics  solutions  are  used  as  a  reference  to  determine  the  accuracy  of  multi-scale  simulations. 
Three  different  approaches  are  used  to  solve  this  problem.  They  are:  (a)  a  macroscopic  model  using  the 
continuum  damage  model  for  constitutive  behavior,  (b)  a  detailed  micromechanical  VCFEM  analysis,  and 
(c)  a  multi-level  model  for  multi-scale  analysis.  The  starting  mesh  in  the  multi-level  model  of  the  bonded 
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joint  consists  of  a  uniform  grid  of  470  QUAD4  elements  for  macroscopic  analysis  as  shown  in  figure  6.11.  The 
constitutive  relation  for  each  element  is  a  fourth  order  anisotropic  CDM  model  that  has  been  developed  for 
this  unit  cell  with  interfacial  cohesive  zone  in  [84],  Figure  6.12(a)  shows  the  contour  of  degrading  dissipation 
energy  Wj,  at  the  final  stage  of  loading  by  a  pure  CDM  based  macroscopic  analysis.  Damage  initiates  near 
the  bottom  left  corner  A  of  the  adhesive  joint  and  propagates  downwards  to  span  the  entire  region  on  the 
left  of  point  A.  Level-0  — >  level- 1  transition  in  the  multi-level  analysis  is  performed  using  equation  (6.22)  and 
level-1  —*  level-2  transition  uses  equation  (6.28)  with  factors  C2  =  0.5  and  C3  =  0.1.  The  gradient  of  the 
energy  +  (75^)^  at  the  final  loading  stage,  used  in  equation  (6.22),  is  shown  in  figure  6.12(b).  The 

corresponding  evolution  of  various  levels  in  the  multi-scale  model  is  depicted  in  figure  6.13  at  two  different 
loading  stages.  There  are  7  level- 1  elements  at  87%  of  the  final  loading.  At  the  final  load  increment,  the 
multi-level  mesh  consists  of  446  level-0  elements,  0  level-1  elements,  14  level-2  elements  and  10  transition 
elements.  All  level-2  elements  emerge  in  critical  the  regions  where  both  the  gradient  and  intensity  of  W<i 
are  high  in  the  macroscopic  analysis.  Figure  6.14(a,b)  depict  the  contours  of  microscopic  stress  an  and  the 
regions  of  debonding  obtained  by  pure  micromechanical  and  the  multi-level  models.  The  results  of  the  multi¬ 
level  model  are  in  excellent  agreement  with  the  micromechanical  analysis,  both  with  respect  to  debonding 
regions  and  evolving  variables.  The  maximum  error  in  an  is  around  1%.  The  excellent  agreement  is  further 
corroborated  in  the  plot  of  an  along  the  vertical  line  through  the  microstructure  in  figure  ??.  Figure  ??(a,b) 
plot  the  macroscopic  (averaged)  En  —  en  curve  obtained  from  (a)  macroscopic  CDM-based  analysis,  (b) 
micromechanical  analysis  and  (c)  multi-scale  analysis  with  the  multi-level  model  at  two  different  locations, 
PI  and  P2  shown  in  figure  6.11(b).  At  P2,  where  the  damage  and  its  gradient  are  low,  solutions  by  the  CDM 
model  and  micromechanics  are  in  relatively  good  agreement.  At  this  point,  the  multi-scale  model  uses  the 
CDM  constitutive  law.  However,  the  CDM  results  are  quite  different  from  the  other  two  at  PI,  a  hotspot 
where  the  damage  and  its  gradient  are  high.  It  is  assuring  to  note  that  the  multi-level  model  matches  the 
micromechanics  results  quite  well  at  this  point. 

The  computational  efficiency  of  the  multi-level  model  is  examined  by  a  comparison  of  the  CPU  time 
on  a  IA32  computer  cluster  for  the  different  models.  The  computations  are  carried  out  in  a  serial  manner 
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using  a  single  processor.  The  results  are  tabulated  in  table  6.1.  Although  the  macroscopic  CDM  analysis  is 
faster,  it  can  lead  to  significant  errors.  The  complete  level-1  solution  is  even  slower  than  the  micromechanics 
solution,  since  it  solves  the  RVE  problem  in  every  element.  Accurate  analysis  with  the  multi-level  model  is  at 
least  7  times  faster  than  the  complete  micromechanics  and  level-1  solutions  for  this  problem.  The  efficiency 
increases  rapidly  with  increasing  number  of  fibers  in  the  analysis. 


Model 

Level-0 

Level- 1 

Micromechanics  (Level-2) 

Multi-scale 

Time  in  seconds 

71 

300330 

300310 

42260 

Table  6.1:  CPU  time  on  a  IA32  cluster  to  solve  the  double  lap  joint  model  by  various  methods. 

Multi-level  analysis  for  model  with  192000  fibers 

This  is  a  more  realistic  model  of  the  composite  joint  with  a  large  number  of  fibers,  to  realize  the  potential 
of  the  multi-level  model.  The  top  ply  consists  of  160  rows  of  fiber,  while  the  bottom  row  consists  of  80  rows 
resulting  in  a  total  of  192000  fibers.  The  geometric  and  material  parameters  are  the  same  as  in  the  previous 
example,  except  for  the  special  cases  mentioned.  A  pure  micromechanical  analysis  is  not  conducted  due 
to  the  large  number  of  fibers.  The  problem  is  analyzed  by  (a)  a  macroscopic  model  by  CDM  and  (b)  the 
multi-level  model.  The  multi-level  analysis  activates  all  three  types  of  adaptation: 

•  Refinement  of  level-0  elements  by  h-adaptation  in  accordance  with  equation  (6.21),  for  C i  =  0.7. 

•  Transition  from  level-0  to  level-1  elements  in  accordance  with  equation  (6.22),  with  C2  =  0.5. 

•  Transition  from  level-1  to  level-2  elements  in  accordance  with  equation  (6.28),  with  C3  =  0.1. 

The  effects  of  variation  of  cohesive  zone  parameters  and  the  effect  of  volume  fraction  are  studied.  The  unit 
cells  considered  in  this  example  have  2  volume  fractions:  (i)  Vf  =  20%.  and  (ii)  Vf  =  40%.  Three  different 
cases  with  different  parameters  in  the  bilinear  cohesive  law  are  considered. 

•  Cl:  Same  cohesive  parameters  as  in  section  6.5.2. 

•  C2:  amax  and  Se  are  the  same  as  in  section  6.5.2.  However,  6C  is  4  times  that  in  case  Cl.  This  reflects 
the  same  cohesive  energy  with  a  smaller  ascending  slope. 


166 


•  C3:  Umax  is  reduced  by  half  and  Se  is  doubled.  Hence  the  cohesive  energy  is  the  same  as  Cl  with  a 
smaller  peak  stress.  Also  6C  is  the  same  as  that  in  Cl. 

The  starting  mesh  has  470  level-0  elements.  For  Vf  =  40%  and  case  Cl,  the  final  mesh  has  1688  level-0 
elements,  24  level-2  elements  and  33  transition  elements  as  shown  in  figure  6.15(a).  Figure  6.15(b)  illustrates 
the  corresponding  microscopic  stress  distribution  and  debonding  in  the  level-2  regions  near  the  hotspot  at 
A.  The  macroscopic  (averaged)  stress-strain  plots  are  shown  for  two  points  in  the  composite  joint:  (a)  near 
the  critical  point  A  and  (b)  at  a  non  critical  point  B  are  shown  in  figure  6.16.  The  predictions  of  the  CDM 
model  agree  with  the  multi-level  model  at  the  point  B.  However,  the  stress  predictions  by  the  CDM  model 
are  considerably  higher  than  those  by  the  multi-level  model  at  A,  where  damage  is  very  localized  and  the 
periodicity  condition  imposed  by  the  CDM  model  is  unrealistic. 

The  effect  of  Vf  on  the  damage  evolution  near  the  corner  PI  is  seen  in  figure  6.17  for  the  case  Cl. 
A  significantly  higher  Wd  is  observed  for  the  higher  volume  fraction,  which  increases  with  evolving  strain. 
Figure  6.18  shows  the  distribution  of  Wd  at  the  end  of  the  analysis  for  the  different  cohesive  parameters. 
Intense  damage  localization  takes  place  near  the  junction  A  in  the  bond  (see  figure  6.11(b).  Damage  starts 
from  this  location  and  propagates  down  and  left  towards  the  edge  of  the  applied  loading.  Damage  localization 
is  the  strongest  for  the  case  Cl,  and  propagates  almost  vertically  down  in  a  narrow  zone.  It  is  in  these  regions, 
that  scale  transition  to  level-2  occurs.  The  damage  distribution  in  the  remaining  parts  of  the  composite  joint 
is  rather  low  and  uniform.  Moving  the  peak  stress  in  case  C2  with  a  lower  traction-displacement  slope  results 
in  a  more  diffused  damage  region  and  the  damage  seem  to  spread  more  in  the  region  to  the  left  of  point 
A.  The  damage  localization  reduces  for  the  case  C3  with  lower  peak  stress  and  the  damage  is  more  evenly 
distributed.  For  Vj  =  20%,  the  damaged  regions  are  less  localized. 

6.6  Conclusions 

An  adaptive  concurrent  multi-level  computational  model  is  developed  in  this  paper  for  multi-scale  analysis 
and  prediction  of  damage  in  fiber  reinforced  composite  materials.  Microstructural  damage  is  manifested  by 
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fiber-matrix  interfacial  debonding  in  this  paper.  Microstructural  damage  mechanisms  leading  to  complete 
failure  are  more  complex  than  the  singular  mode  of  damage  considered  in  this  paper.  The  authors  are 
currently  working  towards  this  goal,  where  interfacial  debonds  bifurcate  into  the  matrix  and  eventually 
coalesce  to  cause  a  continuous  fracture  path.  A  step  forward  in  this  direction  can  be  seen  in  a  recent  paper 
on  the  growth  and  coalescence  of  multiple  cohesive  cracks  [59].  However,  the  intent  of  the  present  paper  is  to 
create  a  framework  for  the  multi-scale  coupling  so  that  more  complex  damage  mechanisms  may  eventually 
be  incorporated.  Hence  interfacial  debonding  is  deemed  sufficient  for  this  purpose. 

The  multi-level  model  invokes  two-way  coupling  of  scales,  viz.  a  bottom-up  coupling  with  homogenization 
at  lower  scales  to  introduce  reduced  order  continuum  models  and  a  top-down  coupling  at  critical  hotspots 
to  transcend  scales  for  following  the  microstructural  damage  evolution.  The  bottom-up  coupling  results  in  a 
continuum  damage  mechanics  (CDM)  model  developed  in  a  preceding  paper  [84].  Three  levels  of  hierarchy, 
with  different  resolutions,  evolve  in  this  model  with  adaptation.  Adaptive  capabilities  enable  effective  domain 
decomposition  in  the  evolving  problem  with  damage,  keeping  a  balance  between  computational  efficiency  and 
accuracy.  Macroscopic  analysis  is  done  with  the  CDM  model  of  [84]  for  high  efficiency.  Pure  micromechanical 
analysis  is  computationally  exhaustive  and  the  adaptive  methodology  optimally  reduces  this  region  to  a 
minimum.  The  Voronoi  cell  finite  element  model  [37,  58]  is  effectively  utilized  for  efficient  micromechanical 
analysis  of  extended  microstructural  regions.  The  numerical  examples  establish  the  accuracy  and  efficiency 
aspects  of  the  model,  as  well  as  demonstrate  its  capability  in  handling  problems  involving  damage  in  large 
composite  domains.  Overall  this  work  lays  an  effective  foundation  for  solving  multi-scale  problems  involving 
localization,  damage  and  crack  evolution  that  may  be  impossible  to  achieve  using  any  single  scale  model. 

6.7  APPENDIX 

6.7.1  Microscopic  and  Macroscopic  Equations  in  Computational  Subdomain 
Level-1  (On ) 

Any  function  /  in  the  RVE  is  assumed  to  be  T-periodic,  i.e.  /(x,  y)  =  /(x,  y  +  kY)  Vk  =  1, 2,  ■  •  • .  Peri¬ 
odicity  conditions  are  used  on  the  RVE  boundary  to  decouple  the  set  of  equations  at  different  levels  as: 
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Microscopic  equations 


fij(x,y)  = 


°u(x>y)  = 

dcrij(x,y) 

dyj 


e  M  )  VUl(y)  |  du>W) 

M  J  +  2(  dyj  +  %  j 


+e/ti1(x) 


dui(  y) 

%  J 


r,  r  X  ,  -1/  N^Wn(y) 

E'ijklOmkOnl  "b  ^ mn 


(Kinematics) 

(x)  (Constitutive) 


0  (Equilibrium) 


(6.33) 


Macroscopic  equations 


Sij(x)  =  | Y\  Jy  Eijki(SkmSin  +  ekJn^)dYemn  =  E^mnemn{^)  (Constitutive) 

~~  +  fi  =  0  (Equilibrium)  (6.34) 


In  the  above  equations  tq  is  a  K-periodic  displacement  function  and  Oij(x,y )  is  the  stress  field  in  the  RVE 
respectively,  while  £;_,(x)  and  eq-  are  the  homogenized  stress  and  strain  tensors.  £qq  and  E^kl  correspond 
to  microscopic  and  homogenized  anisotropic  elasticity  tensor  respectively.  The  details  of  the  derivation  of 
equations  (6.33)  and  (6.34)are  discussed  in  [89,  83]. 


6.7.2  The  Voronoi  Cell  Finite  Element  Model 

A  typical  level-2/  transition  element  consisting  of  microstructural  Voronoi  cell  elements  is  shown  in  figure 
6.3.  The  following  assumptions  are  made  in  the  formulation  of  each  VCFEM  element. 

•  Stress  fields  a”1  in  the  matrix  phase  flm  and  ofj  in  the  inclusion  phase  of  each  Voronoi  cell  element 
are  independent  and  equilibrated.  The  stress  interpolations  in  each  phase  are  expressed  as  : 


{a  m}  =  [P m\{f3m}  in  nm  and  {a  c}  =  [PC]{/3C}  in 


(6.35) 


where  the  matrices  [Pm/C]  are  obtained  from  assumed  stress  functions  like  the  Airy’s  stress  function 
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and  {/3  m/c}  are  unknown  coefficients  to  be  solved. 


•  Compatible  displacement  field  u\  are  assumed  on  each  Voronoi  cell  element  boundary  di \t  and  inter¬ 
polated  as: 

{ue}  -  [Le]{qe}  (6.36) 

•  Compatible  displacement  fields  u™  and  are  assumed  on  the  matrix  and  inclusion  sides  of  the  matrix- 
inclusion  interface  dflc  respectively,  and  are  interpolated  as  : 

{um}  =  [Lc]{qm}  on  dfl™  and  {uc}  =  [Lc]{qc}  on  dClcc  (6.37) 


In  an  incremental  formulation,  the  potential  energy  functional  for  each  element  is  expressed  in  terms  of  the 
incremented  stresses  and  displacements  as: 


ne(a£,  AaVJ,  <4,  Acrfj, «?,  A uf,  <\  A<\  Auf)  =  -  [  A Aa%)dQ 

JQm 

-  [  ABc(afjt  Aa&dfl  +  [  (<r£  +  Aa^  +  A^)dQ 


+ 


r  r  /-u^+au 

/  (4+Aa?.)(4  +Aefi)dn-  /  /  -  ucn)ddn 

Juc  JancJu^- 


-u£ 


r  /‘U^+Au^-uf-Au*  p 

-  /  T?d(uT-uct)ddSl-  (U  +  Ati)(u‘  +  Au')dr 

JdClc  Jl 'tm 


(6.38) 


Here  B  =  \SijkiVij(7ki  is  the  complementary  energy  density  and  A B  =  -S^kiAa^Aaki  +  SljkiA(jl]Cki-  The 
strain  fields  e™  and  are  in  the  matrix  and  inclusion  phases  respectively  of  each  Voronoi  element,  t  is  the 
prescribed  traction  on  the  boundary  Ttm.  The  prefix  A  corresponds  to  increments  and  subscripts  n  and  t 
correspond  to  the  normal  and  tangential  directions  at  the  matrix-inclusion  interface.  The  two  terms  on  the 
matrix-inclusion  interface  8Qm/d flc  provide  the  work  done  by  the  interfacial  tractions  Tm  =  T™ nm  +  T™tm 
due  to  interfacial  separation  (um  —  uc).  T™  and  Ttm  are  the  normal  and  tangential  components  that  are 
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described  by  cohesive  laws  in  [37,  58].  Using  divergence  theorem,  the  potential  energy  can  be  written  as: 


S™klAa™AoS  dQ  - 


[  SSuvZAcr*  dQ 
Jnm 


—  f  -SijklAcr^Aakl  dQ  —  f  SijktahAa^j  dQ 
J  Qc  J  Qc 

+  [  (a%  +  A^)ne«  +  Au')ddQ-  [  (a$  +  Aa^)ncAu^  +  Av?)ddQ 

JdQe  J an™ 


+ 


r  r  /■uir+Au”-<-A< 

/  (4  +  A 4 )n[ ?(u?  +  A ut)ddQ  -  /  /  T; d(u™  -  ucn)ddQ 

J dQ.c  Jd Qc  «/ 


/•up-f-AtiJ”  —  tij  —  Auj 


-  /  /  ' 

JdClc  Ju™—  u® 

-  /  (ti  +  Ati)(uf  +  Auf)dr 

j  rim 


(6.39) 


(6.40) 


Here  ne  and  nc  are  the  outward  normal  on  dfle  and  3QC  respectively.  The  integration  over  the  incremental 
displacements  at  the  interface  dQc  is  conducted  by  the  backward  Euler  method.  The  total  potential  energy 
functional  for  each  level-2  or  tr  element  containing  Nvc  Voronoi  cell  elements  as  shown  in  figure  6.3  is 


N„ 


nl2/tr  =  ne 


(6.41) 


Substituting  stress  interpolations  (6.35)  and  displacement  interpolations  (6.36,6.37)  in  equation  (6.40)  and 
setting  variations  with  respect  to  the  stress  coefficients  Af3m ,  A Pc  respectively  to  zero  results  in  the  weak 
form  of  the  element  kinematic  relation. 


/nm  [pm]T[Sm](pm]dn 

[0] 

< 

’  > 

(3m  +  A/3m 

[0] 

/qc[Pc]T[Sc][Pc]rfn 

pc  +  A(3C 
<  > 

/anJPm]>el[Le]^«  -/anc[P"]>‘][Lc]d«9ft  [0] 

[0]  [0]  Jdnc[P^}T[n^]ddQ 


qe  +  Aqe 
qm  +  Aqm  1 
qc  +  Aqc 


or  in  a  condensed  form 


[He]{/3  +  A/3}  =  [Ge]{q  +  Aq} 


(6.42) 
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The  weak  forms  of  the  global  traction  continuity  conditions  are  subsequently  solved  by  setting  the  variation 
of  the  total  energy  function  in  equation  (6.41)  with  respect  to  Aq,  Aqm  and  Aqc  to  zero.  This  results  in 
the  weak  form  of  the  traction  reciprocity  conditions  as: 


/anjLe]T[nT[Pm]ddft  [0] 

Nvc 

E  “  fan'  [Lc]T[nc]T[Pm]ddfl  [0] 

e=l  m 

[0]  /ang[L‘FKnPc]^ 


0  m  +  A0  m 
0C  +  A/3  c 


> 


Nvc 

E 

e=l 


/rtm  [Le]T{t  +  At}dQ 

-  Jan^  (LC]T  ({nc}T™(w„  +  Aun,  ut  +  A ut)  +  {tc}Ttm(u„  +  A u„,  ut  +  A ut))  ddfl 
~  Jan^F  ({nC}T”(un  +  Aun,  ut  +  A ut)  +  (tc}Ttm(un  +  Aun,  ut  +  Aut))  ddQ. 


(6.43) 


or  in  a  condensed  form 


Nvc  Nvc 

£[Ge]r{/3  +  A/3}  =  £{Re} 


e=l 


e=l 


(6.44) 


Substituting  (6.42)  in  (6.44)  yields: 


i»ve  ■*vvc 

^[Ge]T(He]-1[Ge]{q  +  Aq}  =  £{Re}  (6.45) 

e— 1  e=l 


In  an  iterative  solution  of  equation  (6.45),  its  linearized  form  for  the  i — th  iteration  is  given  as: 


e]T[He]-1[Ge])’{  Aq}’  =  £{Re}*'  -  £[Ge]r[He]-x[Ge]{q  +  Aq}x 

e=l  e.—  1  e=l 


(6.46) 


or  in  a  condensed  form 


[K]'{Aq}’  =  {AF"C}' 


(6.47) 


In  order  to  incorporate  this  relation  in  the  linearized  form  of  the  principle  of  virtual  work  of  equation  (6.11), 
it  should  be  noted  that  the  displacement  vector  u;2/tr  on  the  boundary  of  level-2  and  transition  element  is 


172 


a  subset  of  all  the  VCFEM  displacement  fields,  i.e.  uvc  =  ue  U  um  U  uc.  Consequently,  the  displacement 
field  uuc  can  be  separated  into  two  categories,  viz.  (i)  u,2/tr  on  nodal  points  at  the  boundary  of  the  level-2 
or  transition  elements  shown  is  figure  6.3,  and  (ii)  umt  on  all  the  internal  nodes.  The  stiffness  matrix  and 
the  load  vector  of  the  ensemble  of  all  Voronoi  cell  elements  belonging  to  a  level-2  or  transition  element  can 
therefore  be  partitioned  as 


J£l2/tr,!2/tr  j^L2/tr,int 

t 

Aq12/tr 

t 

r  > 
^pl2/tr 

< 

.  =  < 

j£int,12/tr  j^int.int 

Aqint 

AFint 
«  > 

Static  condensation  of  the  internal  degrees  of  freedom  yields 


^J^l2/tr,l2/trj  _  Jj^l2/tr,intj  jj^int.intj  1  Jj£int,12/trj  j  *  j  ^q12/tr  |* 

=  |AF12/,tr|:  -  [Kl2/tr'int]  [Kint'int]-1  |AFint}1 


(6.48) 


(6.49) 


These  stiffness  matrices  and  load  vectors  are  then  used  in  global  assembly  process  of  equation  (6.11). 
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the  simulation,  for:  (a)  pure 
1-2  region  ( =  0.35),  and 


Transition 


•  Level-2 

(b) 


Figure  6.7:  Computational  mesh  for  the  computational  domain:  (a)  Macroscopic  mesh  with  different  RVE 
in  every  element,  (b)  Multi-level  model  with  the  interface  between  macroscopic  and  microscopic  VCFE 
elements. 


SERVE  boundary 


Level-1  element  boundary 


Figure  6.8:  (a)  A  periodic  microstructure  containing  the  tessellated  RVE  (fibers  in  black),  (b)  placement  of 
the  RVE  in  the  level- 1  element  showing  periodic  nodes  on  the  boundary. 


Figure  6.9:  Comparison  of  microscopic  stress  an  by  different  methods,  plotted  along  a  line  through  the 
middle  of  microstructure 
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Figure  6.10:  Comparison  of  macroscopic  (volume 
at  (a)  macroscopic  element  1,  and  (b)  macroscop 


Micromechanics  (averaged) 
Multi-level  Model  (easel) 
Macroscopic  Level-1 
Multi-level  Model  (casc2) 


0.0009 


0.0012 


veraged)  En  —  cu  curves  by  different  methods  of  analysis 
:  element  9. 


6h 


(a)  (b) 

Figure  6.11:  (a)  Schematic  diagram  of  a  composite  double  lap  joint  showing  dimensions  and  boundary- 
conditions,  (b)  the  level-0  computational  mesh. 
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Figure  6.12:  Contour  plot  of  (a)  degrading  dissipation  energy  Wd  and  (b)  its  gradient  2  +  (f^f)2 

at  the  final  loading  stage. 
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INCREMENT  13/15 
TOTAL  ELEMENT  470 
□  LEVEL  0=  461 
ID  LEVEL  1  =  7 

■  TRANSITION  ELEMENTS  =  5 

■  LEVEL  2  =  4 


(a) 

INCREMENT  15/15 
TOTAL  ELEMENT  470 
□  LEVEL  0=  446 
0  LEVEL  1 =  0 

■  TRANSITION  ELEMENTS  =  10 

■  LEVEL  2=  14 


(b) 

Figure  6.13:  Evolution  of  the  multi-level  computational  model  with  level  transition  (a)  at  87  %  loading,  and 
(b)  at  the  final  loading  stage. 
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Figure  6.14:  Level  2  microscopic  VCFEM  elements  near  the  corner  A  showing  microscopic  stress  distribution 
(GPa)  and  interfacial  debonding  at  the  end  of  the  analysis  by:  (a)  pure  micromechanical  analysis  (b)  multi¬ 
scale  analysis. 
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INCREMENT  IS/ 18 
TOTAL  ELEMENT  1745 
□  LEVEL  0  =  1688 
H  LEVEL  1=  0 

■  TRANSITION  ELEMENTS  =  33 

■  LEVEL  2=  24 


Figure  6.16:  Macroscopic  averaged  stress-strain  (£1 
critical  region  A  and  (b)  non-critical  region  B. 


Damage  Work 


Figure  6.17:  Degrading  dissipation  energy  evolution  near  the  corner  A  of  the  double  lap  joint  for  Vf  =  20% 
and  Vf  =  40%,  and  case  Cl. 
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